{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.fftpack import fft,ifft\n",
    "import pandas as pd\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# N=20取bin的函数\n",
    "\n",
    "def databin_20(lst):\n",
    "    result = [[]]    \n",
    "    length = len(lst)\n",
    "    n = 0\n",
    "    for i in range(length):\n",
    "        result[-1].append(math.log(lst[i],10))\n",
    "        n = n+1\n",
    "        if n == 20:\n",
    "            n = 0\n",
    "            result.append([])\n",
    "    output=[]\n",
    "    for j in range(len(result)):\n",
    "        output.append(np.mean(result[j])) \n",
    "    return output\n",
    "\n",
    "def databin_20_std(lst):\n",
    "    result = [[]]\n",
    "    length = len(lst)\n",
    "    n = 0\n",
    "    for i in range(length):\n",
    "        result[-1].append(math.log(lst[i],10))\n",
    "        n = n+1\n",
    "        if n == 20:\n",
    "            n = 0\n",
    "            result.append([])        \n",
    "    output=[]\n",
    "    for j in range(len(result)):\n",
    "        output.append(np.std(result[j]))\n",
    "    return output\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\npnum = np.arange(len(counts))\\nt = [i*dt for i in pnum]\\n\\n\\nplt.figure(figsize=(10,8))\\nplt.plot(t,counts,\\'b\\')\\nplt.xlabel(\"time\")\\nplt.ylabel(\"counts\")\\nplt.title(\"lightcurve\")\\nplt.show() \\n'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#光变曲线\n",
    "\n",
    "# 提取数据长度\n",
    "data = pd.read_csv(\"0.3_10_final.csv\")  \n",
    "counts_data = data['RATE']\n",
    "dt=100\n",
    "mean_x_data = np.mean(counts_data)\n",
    "\n",
    "N=len(counts_data)\n",
    "N_randomlc=len(counts_data)\n",
    "\n",
    "omega = []\n",
    "POW = []\n",
    "DFT = []\n",
    "fr = []\n",
    "fi = []\n",
    "f1 = []\n",
    "f2 = []\n",
    "p = []\n",
    "f = []\n",
    "\n",
    "\n",
    "counts_list=[None for i in range(500)]\n",
    "\n",
    "# 500条光变曲线\n",
    "for a in range(500):\n",
    "    f_b=1.7E-4\n",
    "    alpha_H=3.8\n",
    "    alpha_L=1.0\n",
    "    for j in range(1,int(N_randomlc)+1):\n",
    "        omega.append(j/(N_randomlc*dt))\n",
    "        POW.append(((omega[-1]**(-alpha_L))/(1+(omega[-1]/f_b)**(alpha_H-alpha_L)))*0.005)\n",
    "        DFT.append(complex(np.sqrt(POW[-1]),np.sqrt(POW[-1])))\n",
    "        s1=np.random.normal(loc=0.0, scale=1.0, size=None)\n",
    "        s2=np.random.normal(loc=0.0, scale=1.0, size=None)\n",
    "        fr.append((DFT[-1].real)*s1)\n",
    "        fi.append((DFT[-1].imag)*s2)\n",
    "        f1.append(complex(fr[-1],fi[-1]))\n",
    "    counts = ifft(f1)\n",
    "    counts_list[a]=counts\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "pnum = np.arange(len(counts))\n",
    "t = [i*dt for i in pnum]\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(t,counts,'b')\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"counts\")\n",
    "plt.title(\"lightcurve\")\n",
    "plt.show() \n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# ????    \n",
    "# 可以直接在光变曲线上乘上数据光变曲线的平均值，相当于加上泊松噪声    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAHwCAYAAADq0mgNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XmcneP9//HXZ5IIg1iKWmfGXrvKIMQelFK0SDB2mhC+StGW0dJWWqXVKu2X2JexJIraiqZFW7UlvraW/pRmYt/3QSq5fn/cJzIZM5kzyZy5zznzej4e53HOvZwzn7lzJG/Xdd3XFSklJEmSlI+avAuQJEnqzwxjkiRJOTKMSZIk5cgwJkmSlCPDmCRJUo4MY5IkSTkyjEkqSxFxWUScXuKf8Y+I2KaUP0OSumMYk9RvpZTWSSndAxARp0XEVTmXVLSImBoR2+ddh6T5ZxiT1Kci4989JRYRA/OuQVJx/AtR6gcKrSgnRcQ/I+LtiLg0IhZsd/ybEfHviHgrIm6OiOUL+38YEecWXg+KiA8j4szC9kIR8XFELFHYHhYRf4+IdyLisfbdfxFxT0SMi4j7gDZglU5q/HJEPBIR70fEdcCCHY7vGhGPFj7/7xGxfoff74SIeDwi3o2I62b9fhGxVETcWnjfWxHx11lhcFbrUkTsBJwMjIqIDwr17x0RUzrUcHxE3NTFNT44Ip4r1P+fiGhqt/++iDi3UNvTETGi3fsWi4iLI+LliHgxIk6PiAEd/myeKnzuPyNio4i4EqgDbinU+52IaIiIFBGHRcQ04M8RsU1EvNDJd2H7wuvTImJiRFxV+PwnImKNwnfltYh4PiJ27Oz3ldR7DGNS/9EEfAVYFVgDOAUgIrYDfgqMBJYDWoFrC++5F9im8Hpj4BVg68L2ZsC/UkpvR8QKwG3A6cCSwAnA7yJi6XY//wBgNLBo4Wd8JiIWAG4Criy8fyKwZ7vjGwGXAGOALwAXADdHxOB2HzMS2AlYGVgfOLiw/3jgBWBp4ItkoWuOdeBSSncAPwGuSyktklLaALgZWDki1mp36v6FGucQEQsDvwZ2TiktCmwOPNrulE2B54ClgFOBGyJiycKxy4FPgdWALwM7AocXPndv4DTgQGAIsBvwZkrpAGAa8LVCvWe2+1lbA2uR/VkX42uF32kJ4P+AO8n+bVgB+BHZtZZUQoYxqf84L6X0fErpLWAcsG9hfxNwSUrpkZTSJ8BJwGYR0QDcD6weEV8AtgIuBlaIiEXI/tG/t/AZ+wO3p5RuTynNTCn9EZgMfLXdz78spfSPlNKnKaX/dqhtGDAI+FVK6b8ppeuBh9sd/yZwQUrpwZTSjJTS5cAnhffN8uuU0kuF3+8WYMPC/v+Shcz6wmf/NRWxKG/hWlxX+N2IiHWABuDWLt4yE1g3IhZKKb2cUvpHu2OvtfvdrgP+BewSEV8EdgaOTSl9mFJ6DfglsE/hfYcDZ6aUHk6Zf6eU5giynTit8Fkfdfc7Fvw1pXRnSulTshC8NHBG4c/oWqAhIhYv8rMkzQPDmNR/PN/udSuwfOH18rRrqUopfQC8CaxQ+Ad9Mlnw2oosfP0dGM6cYawe2LvQFfhORLwDbEEWgjr7+R0tD7zYISS1Dx31wPEdPn+ldr8DZK12s7QBixRenwX8G7ir0I34vbnU0dHlwH4REWQtexMKIW0OKaUPgVHAEcDLEXFbRHyp3Smd/W7LF36vQYX3zPq9LgCWKZy3EvBsD+qFuV/nzrza7vVHwBsppRnttmH2tZRUAoYxqf9Yqd3rOuClwuuXyEIB8FmX2xeAFwu77gW2I+tCe7iw/RVgE+AvhXOeB65MKS3e7rFwSumMdj9zbq1RL5O1uEWHGmd5HhjX4fNrU0rXdPdLp5TeTykdn1JahaxL7tvtx2zNrb6U0gPAdGBLYD866aJsd+6dKaUdyALo08CF7Q539ru9VPi9PgGWavd7DUkprdPu9161qx9ZxP4PgdpZG4WxaEt/7h2ScmUYk/qPoyJixcJYpZPJuuAArgYOiYgNC2OwfgI8mFKaWjh+L9mYpX+mlKYD95B1n/0npfR64ZyrgK9FxFciYkBELFgYPL5ikbXdTzZu6piIGBgR3yALe7NcCBwREZtGZuGI2CUiFu3ugyMb+L9aIQy9B8woPDp6laxLruPfi1cA5wGfppT+1sXP+GJE7FYIsp8AH3T4GcsUfrdBhXFga5F1674M3AX8IiKGRERNRKwaEbPG5V0EnBARQwu/92oRMSs4v0onN0J08P+ABQvXahDZOMHB3bxHUh8zjEn9x9Vk//A/V3icDpBS+hPwfeB3ZC1UqzJ7zBJk3ZILMbsV7J/Ax+22SSk9D+xOFvJeJ2vROZEi/44phLxvkA26f5usy++Gdscnk40bO69w/N/MHqDfndWBSWQB6X7gt7PmFutgYuH5zYh4pN3+K4F1mUurGNnveTxZa9dbZF24Y9sdf7BQxxtk4/X2Sim9WTh2ILAA2XV9G7ieQvduSmli4fyrgffJbnKYNfD/p8Aphe7NEzorKqX0bqGOi8haOj8ku5lBUhmJIsaxSqpwETEVODylNCnvWipNRCxENgB/o5TSM/Pw/oPJrv0WvV2bpOpgy5gkzd2RwMPzEsQkqRjO0CxJXSi0KAawR86lSKpidlNKkiTlyG5KSZKkHOUSxiLirML6bI9HxI3O7ixJkvqrXLopCwvP/jml9GlE/AwgpfTd7t631FJLpYaGhlKXJ0mSNN+mTJnyRkqp24mWcxnAn1K6q93mA8BexbyvoaGByZMnl6YoSZKkXhQR3a0lC5THmLFDgT/kXYQkSVIeStYyFhGTgGU7OdScUvp94ZxmsiVQWubyOaOB0QB1dXVdnSZJklSRShbGUkrbz+14RBwE7AqMSHMZuJZSGg+MB2hsbHQeDkmSVFVyGTMWETsB3wW2Tim15VGDJElSOchrzNh5wKLAHyPi0Yg4P6c6JEmScpXX3ZSr5fFzJUmSyk053E0pSZLUbxnGJEmScmQYkyRJypFhTJIkKUeGMUmSpBwZxiRJknJkGJMkScqRYUySJClHhjFJkqQcGcYkSZJyZBiTJEnKkWFMqiQtLdDQADU12XNLS94VSZLmUy4LhUuaBy0tMHo0tLVl262t2TZAU1N+dUmS5ostY1KlaG6eHcRmaWvL9kuSKpZhTKoU06b1bL8kqSIYxqRKUVfXs/2SpIpgGJMqxbhxUFs7577a2my/JKliGcakStHUBOPHQ309RGTP48c7eF+SKpx3U0qVpKnJ8CVJVcaWMUmSpBzZMtbO1Klw8slz7ovo/HUlH5v1qKmZc7sv9tXUwIAB8/4YOHDuxxZYYM7H4MHZ84ABSJJUlgxj7Sy2GOy66+ztlDp/XcnHOj5mzpz3fTNm9Py9M2dm75vXx6efzv3Yf/8L06fDJ59kz7NeR3w+qHUW2ma9Xnjh7FFbO+dzV69nPS+2GCy+OAwahCRJRTGMtbPEErDffnlXoVKYMWN2OGsf0jrumz4dPvoom0u1rQ0+/DB7tLXBu+/Cyy/P3m5/bNbrd9+Fd96BBRfMvk9LLJGFs66el1oKllsOll0Wllkma92TJPUv/tWvfmHAAFhooexRainBBx/A229nwayz5+eey57feCMLeK+8kr1eYoksmM0KaO0f9fWw8srZ645d0ZKkymUYk3pZBCy6aPboyXysM2ZkgeyVV7LHrJDW2goPPJA9P/dcFvQaGrJgtvLKsMoqs1+vthosskjJfjVJUgkYxqQyMWAAfPGL2WODDbo+74MPsptNnnsO/vOf7HHPPdnzs89m71933dmPddaBL30p6zqVJJUfw5hUYRZZZHbQ6mjGjCyk/eMf8OSTcMst8NOfZiGtvh7WWw823hiGDYOhQ7ObDiRJ+YrU8fa7MtbY2JgmT56cdxlSxZk+Hf7f/4PHH4eHHsq6PZ94AtZYIwtmw4bBpptm2zXOPihJvSIipqSUGrs9zzAm9U+ffAKPPpoFs1mP99+HESNghx2yR3193lVKUuUyjEnqseefh0mTZj8WWwy23z4LZtttl21LkopTbBizQ0LSZ1ZaCQ45BFpasrs5J06EVVeFCy7Iju28M1x4Ibz2Wt6VSlL1MIxJ6lRNTXZX5/HHwx13wIsvZkHtT3/KxpZtsw38+tdZa5okad4ZxiQVZdFFYeRIuPbabP6z44+HRx6BDTfMujCvuCJbhUCS1DOGMUk9tuCC8LWvwWWXwUsvwdixMGECrLgiHHYY/O1vn18nVZLUOcOYpPkyeDDstRfceiv885/ZBLOjR8Oaa8KvfgXvvZd3hZJU3gxjknrNcsvBiSdmk85edlk2XUZDAxxzDDzzTN7VSVJ5MoxJ6nURsPnm2fiyxx/PVg0YPhx23RXuvtsuTElqzzAmqaRWXBF+8pNsofM99oAxY2DLLeGuuwxlkgSGMUl9ZKGF4PDDs3FlY8fCscdmyzDdequhTFL/ZhiT1KcGDoT99ssWMj/hBGhuhk02gT//Oe/KJCkfhjFJuaipgb33hv/7vyyUffOb8NWvZmPMJKk/MYxJylVNDYwaBU89BTvtlK2DefDB8MILeVcmSX3DMCapLCywwOwpMFZYIZvZ/4wzYPr0vCuTpNIyjEkqK0OGwLhx8OCDcN99sN562Z2XklStDGOSytKqq8Itt8AvfgFHHgl77pktVi5J1cYwJqms7bprNqP/uutmXZcXXuhUGJKqi2FMUtlbcEH44Q+z6S8uvBBGjIBnn827KknqHYYxSRVjvfXg/vthl11g003hnHNg5sy8q5Kk+WMYk1RRBgyA44/PQtm118JXvuI0GJIqm2FMUkVafXX4619hq61go43guuvyrkiS5o1hTFLFGjgQvv99uO02OPVUOPBA+OCDvKuSpJ4xjEmqeBtvDI88AoMGQWMjPPZY3hVJUvEMY5KqQm0tXHwxnHIKbL89XHCBU2BIqgyGMUlVZf/94W9/g9/+Fpqa4MMP865IkubOMCap6qy5JjzwQNZtudlm8O9/512RJHXNMCapKi20EFx2GRxxBGy+Odx6a94VSVLnDGOSqlYEjB0Lv/99FsrOOMNxZJLKj2FMUtXbbDN48EGYOBEOPRSmT8+7IkmazTAmqV9YYQX4y1/gnXdghx3gjTfyrkiSMoYxSf3GwgvD734Hw4Zlj6efzrsiSTKMSepnamrgZz+Dk0+GrbeGP/0p74ok9XeGMUn90qGHZutZNjXB+PF5VyOpPxuYdwGSlJdttskWG99lF3j+efjRj7I7MCWpL9kyJqlfW331bMb+P/whm/5ixoy8K5LU3xjGJPV7yywDd98Nzz4LI0fCxx/nXZGk/sQwJknAoovCbbfBgAHw1a/Ce+/lXZGk/sIwJkkFgwfDNdfAWmtl48lefTXviiT1B4YxSWpnwAA47zzYfXfYYgv4z3/yrkhStfNuSknqIAJOPRWWXjqbi2zSJFhjjbyrklStDGOS1IWxY2GhhWC77eCuu2DttfOuSFI1MoxJ0lwccggMGgQjRsCdd8L66+ddkaRqYxiTpG7svz8ssADsuGN2x+XQoXlXJKmaGMYkqQgjR2aBbOed4eabs4XGJak3GMYkqUh77JF1We62G9x4IwwfnndFkqqBU1tIUg/ssgtcdRV8/evwwAN5VyOpGhjGJKmHdtwRLrssayGbPBloaYGGBqipyZ5bWvItUFJFsZtSkubBV78KF10Eu4z4iDum/4Yvf9yaHWhthdGjs9dNTfkVKKli2DImSfNot93gtwscx84f38DjrDf7QFsbNDfnV5ikipJLy1hE/BjYHZgJvAYcnFJ6KY9aJGl+7PnmeD7lbXbiDiaxPWvzVHZg2rR8C5NUMfJqGTsrpbR+SmlD4FbgBznVIUnzp66OUUzgLE5kB/7IM6z22X5JKkYuYSyl9F67zYWBlEcdkjTfxo2D2lqauJrTOI0duYsXFlwt2y9JRchtAH9EjAMOBN4Fts2rDkmaL7MG6Tc3881pF/P2YquwY+1k/vKVxVgq38okVYhIqTSNUhExCVi2k0PNKaXftzvvJGDBlNKpXXzOaGA0QF1d3dDW1tZSlCtJveakk2DSJPjTn2DIkLyrkZSXiJiSUmrs9rxShbFiRUQ9cFtKad3uzm1sbEyTJ0/ug6okad6lBGPHwtNPw+23w0IL5V2RpDwUG8ZyGTMWEau329wNeDqPOiSpFCLgvPNg2WVh1Cj473/zrkhSOcvrbsozIuLJiHgc2BH4Vk51SFJJDBgAV1wBM2bAoYfCzJl5VySpXOV1N+WeKaV1C9NbfC2l9GIedUhSKQ0aBBMnZpPyH3NM1n0pSR05A78klVBtLdxyC/z973DaaXlXI6kcuTalJJXYYovBHXfA5pvD8svDmDF5VySpnBjGJKkPLLMM3HknbLllNrB/993zrkhSuTCMSVIfWXVVuPlm2HlnWHrprKVMkhwzJkl9qLERrrwSvvGNbB4ySTKMSVIf22knOOOM7Pmll/KuRlLe7KaUpBwcfHAWxL76Vbj33myQv6T+yZYxScrJSSfB8OFZl+Unn+RdjaS8GMYkKScR8OtfZ61ihx7qpLBSf2UYk6QcDRgALS3w7LNOCiv1V44Zk6ScLbRQNuXFsGGwyipw0EF5VySpLxnGJKkMLLMM3HYbbLMN1NXBttvmXZGkvmI3pSSVibXWgmuvhX32gaeeyrsaSX3FMCZJZWTbbeHMM2GXXeC11/KuRlJfMIxJUpk56CDYf/9s/cqPPsq7GkmlZhiTpDL0wx9mg/kPPBBmzsy7GkmlZBiTpDIUAZdcAq+8AiefnHc1kkrJMCZJZWrwYLjpJrjhBrjwwryrkVQqTm0hSWXsC1+A22+HLbbIui1HjMi7Ikm9zZYxSSpzq62WTXmx337wzDN5VyOptxnGJKkCbLMNnH46fO1r8M47eVcjqTcZxiSpQnzzm7DTTjByJHz6ad7VSOothjFJqiA//znU1MBxx+VdiaTeYhiTpAoycCBcdx1MmgTnn593NZJ6g3dTSlKFWWwxuOWW7A7LNdaA7bbLuyJJ88OWMUmqQKutBtdcA/vu6x2WUqUzjElShdp2W/jRj7zDUqp0hjFJqmBjxsBXvgKjRnmHpVSpDGOSVOF+8YtsLctvfzvvSiTNC8OYJFW4gQOzGfr/+EcYPz7vaiT1lHdTSlIVWHxx+P3vszss11kHhg/PuyJJxbJlTJKqxBprwBVXwN57wwsv5F2NpGIZxiSpiuy0E3zrW/D1r8NHH+VdjaRiGMYkqcp85zvZPGRjxkBKeVcjqTuGMUmqMhFw8cXwxBPwq1/lXY2k7jiAX5KqUG0t3HQTDBsG664LO+yQd0WSumLLmCRVqfr6bMqL/feHZ5/NuxpJXTGMSVIV23pr+MEPYI894IMP8q5GUmcMY5JU5caOhU03hYMOgpkz865GUkeGMUmqchHwm9/Ayy/DuHF5VyOpIwfwS1I/MHgw/O53sMkmsMEGsNtueVckaRZbxiSpn1huObj+ejj8cHjqqbyrkTSLYUyS+ouWFjYd1cCZrx/C7us/x9vjJ+ZdkSQMY5LUP7S0wOjR0NrKwVzGzp/ewn5HLsaMK6/OuzKp3zOMSVJ/0NwMbW2fbf6cE/h45iBOPfrNHIuSBIYxSeofpk2bY3MQn3Ido7jivd256aacapIEGMYkqX+oq/vcrmV4neuX/R9Gj4ann86hJkmAYUyS+odx47IFK9urrWWTn4/kpz+Fr38d3nsvn9Kk/s4wJkn9QVMTjB+fLVgZkT2PHw9NTRx2GGy1FRx8MKSUd6FS/xOpgv7La2xsTJMnT867DEmqOp98kq1jufvucNJJeVcjVYeImJJSauzuPGfglyQxeHA2Iewmm8DQobDjjnlXJPUfdlNKkgBYcUW45ho44AD4z3/yrkbqPwxjkqTPbL01nHwyfOMbc0xLJqmEDGOSpDkccwysvTYccYQD+qW+YBiTJM0hAi68EB57DM47L+9qpOrnAH5J0ufU1sKNN8Jmm8GGG8KWW+ZdkVS9bBmTJHVqlVXg8sthn33gpZfyrkaqXoYxSVKXdtoJxo6FvfaC6dPzrkaqToYxSdJcnXQSLLMMHHts3pVI1ckwJkmaq5qarLvyT3+CSy/Nuxqp+jiAX5LUrcUWywb0b701bLABbLRR3hVJ1cOWMUlSUdZeG37zm2z82Ftv5V2NVD0MY5Kkoo0cmS0mfsABMHNm3tVI1cEwJknqkTPPhPfeg3Hj8q5Eqg6GMUlSjwwaBBMmwPnnw5135l2NVPkMY5KkHltuObjmGjjoIGhtzbsaqbIZxiRJ82SrreCEE7IB/R9/nHc1UuUyjEmS5tnxx0N9vRPCSvPDMCZJmmcRcMklcM892cSwknrOSV8lSfNlyBD43e9gm22yCWE33DDviqTKYsuYJGm+rbMOnHsu7LknvP123tVIlaXbMBYRS/ZFIZKkyrbPPrDrrtkdlk4IKxWvmJaxByNiYkR8NSKi5BVJkirWWWfBm2/CGWfkXYlUOYoJY2sA44EDgH9HxE8iYo3SliVJqkQLLJBNCHveeTBpUt7VSJWh2zCWMn9MKe0LHA4cBDwUEfdGxGYlr1CSVFFWWAFaDpnEATu9xvNRBw0N0NKSd1lS2SpmzNgXIuJbETEZOAH4H2Ap4Hjg6hLXJ0mqNC0tbPur3Tluxs/Zi4l80voyjB5tIJO6UEw35f3AEGCPlNIuKaUbUkqfppQmA+eXtjxJUsVpboa2Nk7kLFbgRb7N2dDWlu2X9DnFzDO2ZkopdXYgpfSzXq5HklTppk0DIIBLOYSNeZiraGL/aXamSJ0ppmXsrohYfNZGRCwREXf2xg+PiBMiIkXEUr3xeZKkMlBX99nLxXiPG/gGx/FLHl92xxyLkspXMWFs6ZTSO7M2UkpvA8vM7w+OiJWAHYBp8/tZkqQyMm4c1NZ+trku/+BXC3yXPWdO4J135vI+qZ8qJozNiIjP/jcnIuqBTrste+iXwHd66bMkSeWiqQnGj89WEI+A+nqaLhnBV/YawsEHOyGs1FExYawZ+FtEXBkRVwJ/AU6anx8aEbsBL6aUHpufz5EklammJpg6NUteU6dCUxNnnw2vvppNDNtjLS3ZFBk1NU6VoarT7QD+lNIdEbERMIxsPOZxKaU3untfREwClu3kUDNwMlDU4IGIGA2MBqhrNw5BklRZZk0Iu8kmMGwYbL11kW9sacmmxmhry7ZbW7NtyEKfVOGiixsl5zwpYgWgnnbhLaX0l3n6gRHrAX8CCv9VsSLwErBJSumVub23sbExTZ48eV5+rCSpTNx1FxxyCEyZAst29r/sHTU0ZAGso/r6rNVNKlMRMSWl1Njded22jEXEz4BRwD+AWT39iay7ssdSSk/Q7gaAiJgKNBbT2iZJqnw77pg1bO2zT7Zk0sDu/iWa1sV9Xl3tlypMMWPG9iCba2yXlNLXCo/dSl2YJKl6nXIKDB4M3/9+ESd3NUTFoSuqEsWEseeAQaUqIKXUYKuYJPUvAwbAVVdlw8FuuaWbkztMlQFk2+PGlaw+qS8VMwN/G/BoRPwJ+GTWzpTSMSWrSpJU9ZZeGq67DvbYAx54AFZeuYsTZw3Sb27Ouibr6rIg5uB9VYluB/BHxEGd7U8pXV6SiubCAfySVH3OOQeuuALuuw8WXDDvaqTeU+wA/mLvplwIqEsp/as3iptXhjFJqj4pwahRsOSScP75eVcj9Z5iw1i3Y8Yi4mvAo8Adhe0NI+Lm+S9RkqRskv6LLoK774Yrr8y7GqnvFTOA/zRgE+AdgJTSo0BXPfuSJPXYkCFw/fXw7W/Dk0/mXY3Ut4oJY5+mlN7tsM/1JCVJvWq99eAXv4C99oL338+7GqnvFBPGnoyI/YABEbF6RJwL/L3EdUmS+qEDD4SttoLDD8/Gkkn9QTFh7H+AdcimtbgGeA84tpRFSZL6r1//Gp55Bs47L+9KpL5RzELhbWSLezeXvhxJUn+34IIwcSJsthlsvHG2qLhUzYpZm/JuOhkjllLariQVSZL6vVVXhQsvzKa8mDIFlloq74qk0ilmBv4T2r1eENgT+LQ05UiSlNl9d/j732H//eH226GmmIE1UgXq9qudUprS7nFfSunbwKZ9UJskqZ8bNw4++ghOPz3vSqTSKaabcsl2mzXAUGDZklUkSVLBwIFw7bUwdGg2hmyHHfKuSOp9xXRTTiEbMxZk3ZP/AQ4rZVGSJM2y3HLQ0gL77QcPPwwrrph3RVLvKuZuSmfblyTlattt4ZhjYORIuPdeGDQo74qk3lNMN+U35nY8pXRD75UjSVLnvvtduO++7Pnss/OuRuo9xXRTHgZsDvy5sL0tcA/wLln3pWFMklRyNTVwxRXZ+LHNN8+WTZKqQTFhLAFrp5ReBoiI5YDfpJQOKWllkiR1sOSS2YSwO+8M668Pa6yRd0XS/Ctm1paGWUGs4FXAr78kKReNjfDjH2ctY21teVcjzb9iwtg9EXFnRBwcEQcBtwF3l7guSZK6NGZM1jI2dqwLiqvyFTPp69HA+cAGwIbA+JTS/5S6MEmSuhIB55+fTXVx8cV5VyPNn2LGjAE8AryfUpoUEbURsWhK6f1SFiZJ0twssghcfz1stVU2qP/LX867ImnedNsyFhHfBK4HLijsWgG4qZRFSZJUjLXWgnPPzcaPvfNO3tVI86aYMWNHAcOB9wBSSs8Ay5SyKEmSirXPPtndlQcf7PgxVaZiwtgnKaXpszYiYiDZdBeSJJWFX/wCXn45e5YqTTFh7N6IOBlYKCJ2ACYCt5S2LEmSijd4MEyYAGedlc3SL1WSYsLY94DXgSeAMcDtwCmlLEqSpJ6qr8/urNxnH3j99byrkYo317spI2IAcHlKaX/gwr4pSZKkebPrrvC3v8H++8Ptt8OAAXlXJHVvri1jKaUZwNIRsUAf1SNJ0nw5/XT46CMYNy7vSqTiFDPP2FTgvoi4Gfhw1s6U0tmlKkqSpHk1cCBce20299jw4TBiRN4VSXNXzJixl4BbC+cu2u4hSVJZWn55uOoqOOAAeOkloKUFGhqgpiZ7bmnJuUJpti5bxiLiypTSAcA7KaVz+rAmSZLm24gRcMQRsO+IV/neDDiFAAAgAElEQVRT65EM/KiwcExrK4wenb1uasqvQKlgbi1jQyOiHjg0IpaIiCXbP/qqQEmS5lVzMwye+i++/9FJcx5oa8sOSmVgbmPGzgfuAFYhW5uyvVTYL0lS2RowAFo+3ouNmMJw7mNXbpt9cNq0/AqT2umyZSyl9OuU0lrAJSmllTs8DGKSpIqwdH0t17Avh3ExrdTNPlBX1/WbpD401wH8EbFfSunIiNinrwqSJKlXjRvHFrX/x4mcxUgmMJ1BUFvr3BcqG93dTblCRIwEVuyLYiRJ6nVNTTB+PMfXXc+yvMKJi14A48c7eF9lo8swFhGnAksCVwNLRsQP+qwqSZJ6U1MT0TqVy97anVuWOoTrBxvEVD7mNmbsh8BbwP7AWymlH/VZVZIklcASS2QLih95JDzzTN7VSJnuuilfSildC7zYF8VIklRqjY3wwx/C3ntnyyZJeesujF0TEU+mlK7pk2okSeoDRx4JX/oSHHNM3pVI3S8UPhN4LCK8/1eSVDUi4MIL4S9/gSuu6OIkl1BSHylmofDlgH9ExEPMuVD4biWrSpKkElt0Ubj+ethuu2xR8XXWaXewpSVbMqmtLdt2CSWVUKSU5n5CxNad7U8p3VuSiuaisbExTZ48ua9/rCSpil16KZx1Fjz0ECyySGFnQ0MWwDqqr4epU/uwOlWyiJiSUmrs7rzuxozNCl1TgUGF1w/z+eWRJEmqSIccAsOGwZgx8Fn7RFdLJbmEkkqg2zAWEd8ErgcuKOxaAbiplEVJktSXzjsPnngimwsW6HqpJJdQUgl0G8aAo4DhwHsAKaVngGVKWZQkSX2pthYmToRTToFHHiFbKqm29vMnuYSSSqCYMPZJSmn6rI2IGAjMfaCZJEkVZs01sxayvfeGd3bJllCivj679bK+3iWUVDLF3E15b0ScDCwUETsAY4FbSluWJEl9b9Qo+Otf4dBD4Xe/ayIMX+oDxbSMfQ94HXgCGAPcDpxSyqIkScrLL36RjdM/55y8K1F/UUzL2DZAS0rpwhLXIklS7gYPzsaPbbpp9thss7wrUrUrpmXsYODRiLg/Is6MiK9FxBIlrkuSpNysvDJcdBHssw+88Ube1ajaFTPP2IEppTWAPYEXgN+QdVtKklS1dtsNRo6EAw+EmTPzrkbVrJh5xvaPiAvI5hrbHjgP2LLUhUmSlLef/ATeew/OOCPvSlTNihkz9ivgWeB84O6U0tSSViRJUpkYNAiuvRY23jgbO7bttnlXpGpUTDflUsChwILAuIh4KCKuLHllkiSVgRVXhMsvh/33h1deybsaVaNiuimHAHVAPdAALAbYey5J6jd23BEOPxz23RdmzMi7GlWbYu6m/BvwNeBxYFRKac2U0kGlLUuSpPLygx/AgAFw6ql5V6Jq0+2YsZTS+gARsSgugyRJ6qcGDICWFhg6FIYPh513zrsiVYtiuinXjYj/A54E/hkRUyJi3dKXJklSefniF+Hqq+GQQ+D55/OuRtWimG7K8cC3U0r1KaU64PjCPkmS+p2ttoLjjsvmIJs+Pe9qVA2KCWMLp5TunrWRUroHWLhkFUmSVOZOPBGWWgq+9728K1E1KCaMPRcR34+IhsLjFOA/pS5MkqRyVVOTTXdxww1w4415V6NKV0wYOxRYGrgBuLHw+pBSFiVJUrlbckmYMAHGjIFnn827GlWyYiZ9fTuldAywLbBVSulbKaW3S1+aJEnlbZNN4Pvfh71HvMXH9WtmTWYNDdltl1KRirmbcuOIeAJ4DHgiIh6LiKGlL02SpPJ39BItrPrCPRw37VhICVpbYfTozgNZS0sW1gxtaqeYbsqLgbEppYaUUgNwFHBpSauSJKlCxCnNXDTjECaxPVezb7azrQ2am+c8saUlC2mtrd2HNvUrxYSx91NKf521kVL6G/B+6UqSJKmCTJvGYrzH9ezFtziHp/jSZ/vn0NychbT2Ogtt6neKCWMPRcQFEbFNRGwdEb8F7omIjSJio1IXKElSWaurA2ADHuennMReXM+H1H62/zMdw1l3+9VvdLscErBh4bnjalybky2PtF2vViRJUiUZNy7rbmxr4zAu5q9syZEDLuTy0xPR/ry6uqxrsqOOoU39TjFrU27bF4VIklSRmpqy5+ZmYto0frviGWya7ueSTxbnsPbntQttn6mtzfarXyumm1KSJM1NUxNMnQozZ7LwtKeYeNfifO978OijHc4ZPx7q6yEiex4/fnaYU79VTDelJEnqgbXWgnPOgb33hilTYMiQwoGmJsOXPseWMUmSSmC//WDECDj88GwmC6kr3baMRcQ3Otn9LvBESum13i9JkqTq8KtfwWabwW9/C0cdlXc1KlfFdFMeBmwG3F3Y3gZ4AFgjIn6UUrqyRLVJklTRFlwQJk6EzTeHTTeFxsa8K1I5KqabciawVkppz5TSnsDawCfApsB3S1mcJEmVbrXVspaxkSPhbVd2VieKCWMNKaVX222/BqyRUnoL+G9pypIkqXrstRfsuisccojjx/R5xYSxv0bErRFxUEQcBNwM/CUiFgbemZcfGhGnRcSLEfFo4fHVefkcSZIqxVlnwUsvZePIpPaKGTN2FPANYAsggMuB36WUEjA/E8L+MqX08/l4vyRJFWPwYJgwATbZJBvUP2xY3hWpXBQzA3+KiL8B08mWP3qoEMQkSVIPNDTARRfBqFHwyCPwhS/kXZHKQbfdlBExEngI2AsYCTwYEXv1ws8+OiIej4hLImKJXvg8SZLK3m67ZZPBHnggzJyZdzUqB9FdI1dEPAbsMGtOsYhYGpiUUtqgm/dNApbt5FAz2dQYb5C1tP0YWC6ldGgXnzMaGA1QV1c3tLWzRVYlSaog//0vbL11Fsy+9728q1GpRMSUlFK3E5oUM2aspsPkrm9SRItaSmn7Ij6biLgQuHUunzMeGA/Q2Nho96gkqeINGgTXXQcbb5zNQbbVVnlXpDwVczflHRFxZ0QcHBEHA7cBt8/PD42I5dptfh14cn4+T5KkSrPSSnDppdmySa+5nk2/VswA/hMjYk9gONndlONTSjfO5889MyI2JOumnAqMmc/PkySp4uy8Mxx0ULZ2+B13wIABeVekPHQ7ZqycNDY2psmTJ+ddhiRJvebTT2H77WG77eAHP8i7GvWm+R4zFhHvk7Vcfe4Q2YwXQ+ajPkmSBAwcCFdfna1bOXw4jBiRd0Xqa12GsZTSon1ZiCRJ/dXyy8OVV8IBB8CUKbDcct2/R9WjmAH8kiSpxEaMgDFjYN99s67LLrW0ZLPH1tRkzy0tfVShSsUwJklSmTjllGzai9NO6+KElhYYPRpaW7MVx1tbs20DWUUzjEmSVCYGDMhy1WWXZXdXfk5zM7S1zbmvrS3br4plGJMkqYwss0w2oP/gg+H55zscnDat8zd1tV8VwTAmSVKZ2WorOPZY2GefbOmkz9TVdf6GrvarIhjGJEkqQ9/5Diy+OJx8crud48ZBbe2cJ9bWZvtVsQxjkiSVoZoauOIKmDABbr65sLOpCcaPh/p6iMiex4/P9qtiOQO/JEll7P77YY894MEHs5ksVDmKnYHfljFJksrYZpvBd78LI0fC9Ol5V6NSMIxJklTmjjsum6X/xBPzrkSlYBiTJKnMRcCll8Itt8D11+ddjXqbYUySpAqwxBLZYP6xY+Hf/867GvUmw5gkSRWisRF+8APYe2/4+OO8q1FvMYxJklROulkI/KijYPXVs3Fkqg6GMUmSykURC4FHwEUXwaRJ2bJJqnyGMUmSykWRC4EPGQITJ8K3vgVPP92H9akkDGOSJJWLHiwEvuGG8JOfZOPHOuY3VRbDmCRJ5aKHC4EffjhssAEcfXQJa1LJGcYkSSoXPVwIPALOPz9bMumyy0pfnkrDMCZJUrmYh4XAF1kkmwj2xBPhySf7sFb1moF5FyBJktppappr+OrMOuvAz3+ejR97+OEsoKly2DImSVIVOOgg2HxzGDMmmxVDlcMwJklSlTj3XHj8L+9w4VIndTlprMqPYUySpCpRe2ML17+xDc1vfZtH0/qdThqr8mMYkySpWjQ3s+bHj/FrjmFvJvIuQzqdNFblxTAmSVK1KEwOuy/Xsj2TOJyLSO32qzwZxiRJqhbtJof9JcfxLKvyG47qejJZlQXDmCRJ1aLdpLEL8gkT2ZsfcSoPH/q/ORemuTGMSZJULTpMGrtq/Qz+95inGHnpzrz9dt7FqStO+ipJUjXpMGnsnsBfAw4+GG66KZvYX+XFljFJkqrcmWfCK6/A2WfnXYk6YxiTJKnKLbAAXLfvTZz5ndf5ewx3MtgyYxiTJKnatbTQ0NzExTMPYR+u4Y3WD5wMtowYxiRJqnbNzdDWxq7cxj5cy4Fcwcy2j5wMtkwYxiRJqnbtJn0dRzPvMYSf8V0ngy0T3k0pSVK1q6vL1qkEBvEp17IPjUxm82WeY+ucS5MtY5IkVb92k8ECrMiLXDb4CJqmX8qrr+ZYlwDDmCRJ1a/DZLDU17PTxXtzyFG1NDXBjBl5F9i/RUop7xqK1tjYmCZPnpx3GZIkVYUZM2D77WGbbeDUU/OupvpExJSUUmN359kyJklSPzVgAFx9NVxwAUyalHc1/ZdhTJKkfmy5P7dw1Yx9OWCHl3lpxU1mzz3W0pJNDltT4ySxJWYYkySpv2ppgdGj2e61axnLb9n3xbP49JtHwtix2aSwra2QUvY8enS234DW6xwzJklSf9XQ8NmUFzMJduYPDGUKPxnwg85H9Udk4WyW2trsxoB2C5NrNseMSZKkuWs36WsNiavYnys5gD/M2KHz8zs24LS1OYt/LzCMSZLUX9XVzbG5NG9wNftxCJfyPCsW9xnO4j/fDGOSJPVXHSaDBdiy9hGOG3Y/o2om8t/2C/VEdP4ZHQKdes4wJklSf9XJZLCMH8+J932dJddfiZOG/Gb2/iOO+Fxwo7Y2C3SaLw7glyRJn/PmmzB0KJxzDuy+e2FnS0s2RmzatKxFbNw4B+/PRbED+A1jkiSpUw88kAWxBx6AlVfOu5rK492UkiRpvgwbBiedBCNHwief5F1N9TKMSZKkLn3rW7DSSnDCCfP5Qc7o3yXDmCRJ6lIEXHIJ3H47TJw4jx9SmOn/czP6G8gAw5gkSerG4ovDhAlw1FHwzDPz8AHNzdkEse05YexnDGOSJKlbQ4fCaadl48c++qiHb+5qYlgnjAUMY5IkqUhHHglrrgnHHtvDN3Y1MawTxgKGMUmSVKSIbI7Yu+/u4XCvTmb6d8LY2QxjkiSpaEOGZAP5jz0WnnqqyDd1MdO/E8ZmnPRVkiT12MUXwy9/CQ8+CAsvnHc15clJXyVJUskceihstBEcfXTelVQ+w5gkSeqxCPjf/4WHHoJLL827mso2MO8CJElSZVp44Wz82NZbQ2MjrLde3hVVJlvGJEnSPFt7bTj7bNh7b3j//byrqUyGMUmSNF8OOAC23BLGjMlWO1LPGMYkSdJ8+/Wv4R//yGasUM8YxiRJ0rxraYGGBhZauIaJb2zLKSd+zP+Nux0aGqCmJnt2QfC5MoxJkqR509ICo0dDayukxBov3cN5bYex9ylr8m7r21mfZWtrdo6BrEuGMUmSNG+am6GtbY5do2ZczVe4g8O4mM+Gj7W1ZeeqU4YxSZI0b6ZN63T32XybqTRwLv/T7bkyjEmSpHlVV9fp7sFMZwIjOZ1TeIiN53pulwpj0frDuDPDmCRJmjfjxkFt7Zz7Bg2CBRZgFf7DBYxhFNfx1kIrZOcWq8NYtGofd2YYkyRJ86apKZvLor4+Wx+pvj5bG+mSS6C+nq/H79lj0T9z8JceIO3XVPzndjIWrZrHnUWqoNnZGhsb0+TJk/MuQ5IkFWn6dNhqK9hrLzjhhCLfVFPT+eyxETBzZq/WV0oRMSWl1NjdebaMSZKkkllgAZgwAc46C+67by4nth8jVtNFPOnpuLMKYRiTJEklVVeX9Vzuuy+8/nonJ3QcIzZjxufPqa3t2bizCmIYkyRJJbfLLrDfftk6lp/raexsjBjAgAGzx6KNH5+NUatChjFJktQnTj8dPvwQzjijw4Gu5iCbOTN7TJ1atUEMcgxjEfE/EfGviPhHRJyZVx2SJKlvDBwI114L554L99zT7kBXY8GqdIxYR7mEsYjYFtgdWD+ltA7w8zzqkCRJfWuFFeDyy7OGrldfLezsbL6yKh4j1lFeLWNHAmeklD4BSCm9llMdkiSpj+24Ixx2WDaGbMYMOp+vrIrHiHWUVxhbA9gyIh6MiHsjYuOc6pAkSTk49dTsxskf/aiwo6kpGxvWD8aIdTSwVB8cEZOAZTs51Fz4uUsAw4CNgQkRsUrqZAbaiBgNjAao6yd9x5IkVbsBA+Dqq2HoUNhiC9hhh7wryk8uM/BHxB1k3ZT3FLafBYallDqbfeQzzsAvSVJ1ufvurLty8uRsPFk1KfcZ+G8CtgOIiDWABYA3cqpFkiTlZNtt4eijswlhP/0072rykVcYuwRYJSKeBK4FDuqsi1KSJFWpdssfnTR+ZWrfeYnvfz/vovJRsjFjc5NSmg7sn8fPliRJOZu1/FFh1v2aaVO5cqHN2Gj8P9lii4XZZZec6+tjzsAvSZL6VifLHy390TSuHXQghx3W9YT8n2m/qHhDQ7ZdwQxjkiSpb3WRtoa/diPHHw+jRsH06e0OtA9fSy0Fhx46e1Hx1tasla2CA5lhTJIk9a25LH90/PGw9NLwve8V9s3q0pwVvt58s0NSI2tla24uacmlZBiTJEl9ay7LH9XUwGWXwQ03wI030mmXZqe67dssX4YxSZLUt7pZ/mjJJeG662DMGHiudUBxn1nBE8MbxiRJUt/rZvmjTTfNGsVGLnAjn7DA3D+rwhcVN4xJkqSydMwxUL/B4hw/8Jw5DwwaBF/4QtUsKm4YkyRJZSkCLvljHXcs2cSEpcbODl+XXgpvvFE1i4obxiRJUtlabDGY8IdFOTp+wzP/6iR8FTPnWJnPS2YYkyRJZW2jjeCHP4S994aPPmp3oOO0F53NOVbMOTmLSloSsrGxMU2ePDnvMiRJUh9LCfbbDxZdNBsiBmStXK2tnz+5vj5rQSv2nBKJiCkppcbuzrNlTJIklb2ILITdey9cdVVhZ1dzi7XfX8w5OTOMSZKkirDoojBxIhx3HDz1FHOdyb/T112dkzPDmCRJqhjrrw9nnJGNH/vw+2d0OZP/Z+Yy23+5MIxJkqSKcuihMHQoHPXXfeY6kz/Q7Wz/5cAB/JIkqeJ8+CFssgmccAIcckje1XSu2AH8A/uiGEmSpN608MLZ+LGtt4bGRlhvvbwrmnd2U0qSpIq09tpw9tnZ+LH338+7mnlnGJMkSRXrgANgiy3giCOyucgqkWFMkiRVtHPPhSefbDcZbIVxzJgkSapoCy2UjR8bPjwb1P/lL+ddUc/YMiZJkireGmvAeedl48fefTfvanrGMCZJkqrCqFGw445w+OGVNX7MMCZJkqrG2WfDc89lrWSVwjFjkiSpaiy4IEyYAJttBsOGwcYb511R92wZkyRJ1aGlBRoaWHX1Gs5nDCN3+YC33867qO4ZxiRJUuVraYHRo6G1FVLiG6+PZ/e3L+fg7Z8v+/FjhjFJklT5mpuhrW2OXWd+ehyvPPkmZ5+dU01FMoxJkqTKN23a53YtwH+ZMH0PzjwT7r8/h5qKZBiTJEmVr66u09319XDxxbDPPvDGG31cU5EMY5IkqfKNGwe1tXPuq62FcePYdddsDrIDDyzP+cec2kKSJFW+pqbsubk567Ksq8sCWmH/uHFw770QkWONXYhUjhGxC42NjWny5Ml5lyFJktStiJiSUmrs7jy7KSVJknJkGJMkScqRYUySJClHhjFJkqQcGcYkSZJyZBiTJEnKkWFMkiQpR4YxSZKkHBnGJEmScmQYkyRJypFhTJIkKUeGMUmSpBwZxiRJknJkGJMkScqRYUySJClHhjFJkqQcGcYkSZJyZBiTJEnKkWFMkiQpR4YxSZKkHBnGJElS9WlpgYYGqKnJnlta8q6oSwPzLkCSJKlXtbTA6NHQ1pZtt7Zm2wBNTfnV1QVbxiRJUnVpbp4dxGZpa8v2lyHDmCRJqi7TpvVsf84MY5IkqbrU1fVsf84MY5IkqbqMGwe1tXPuq63N9pchw5gkSaouTU0wfjzU10NE9jx+fFkO3gfvppQkSdWoqalsw1dHtoxJkiTlyDAmSZKUI8OYJElSjgxjkiRJOTKMSZIk5cgwJkmSlCPDmCRJUo4MY5IkSTkyjEmSJOXIMCZJkpQjw5gkSVKODGOSJEk5MoxJkiTlyDAmSZKUI8OYJElSjgxjkiRJOTKMSZIk5cgwJkmSlKNIKeVdQ9Ei4nWgNecylgLeyLmG/sTr3be83n3Pa963vN59q79f7/qU0tLdnVRRYawcRMTklFJj3nX0F17vvuX17nte877l9e5bXu/i2E0pSZKUI8OYJElSjgxjPTc+7wL6Ga933/J69z2ved/yevctr3cRHDMmSZKUI1vGJEmScmQY60ZEnBYRL0bEo4XHV7s4b6eI+FdE/DsivtfXdVabiDghIlJELNXF8Rnt/kxu7uv6qk0R1/ugiHim8Dior+urFhHx44h4vPC9vSsilu/iPL/fvaQH19zveC+IiLMi4unCNb8xIhbv4rypEfFE4c9lcl/XWW7spuxGRJwGfJBS+vlczhkA/D9gB+AF4GFg35TSP/ukyCoTESsBFwFfAoamlD43R01EfJBSWqTPi6tC3V3viFgSmAw0AgmYUjjv7b6utdJFxJCU0nuF18cAa6eUjujkPL/fvaSYa+53vPdExI7An1NKn0bEzwBSSt/t5LypQGNnf7/3R7aM9Y5NgH+nlJ5LKU0HrgV2z7mmSvZL4Dtkfymq9Lq73l8B/phSeqvwj9MfgZ36qrhqMisUFCyM3/GSK/Ka+x3vJSmlu1JKnxY2HwBWzLOeSmEYK87RhSbXSyJiiU6OrwA83277hcI+9VBE7Aa8mFJ6rJtTF4yIyRHxQETs0Re1VaMir7ff714UEeMi4nmgCfhBF6f5/e5FRVxzv+OlcSjwhy6OJeCuiJgSEaP7sKayNDDvAspBREwClu3kUDPwv8CPyb44PwZ+QfYFm+MjOnmv/8fbhW6u98nAjkV8TF1K6aWIWAX4c0Q8kVJ6tjfrrBa9cL39fvfA3K53Sun3KaVmoDkiTgKOBk7t5Fy/3z3QC9fc73gPdHe9C+c0A58CLV18zPDCd3wZ4I8R8XRK6S+lqbj8GcaAlNL2xZwXERcCt3Zy6AVgpXbbKwIv9UJpVamr6x0R6wErA49FBGTX8ZGI2CSl9EqHz3ip8PxcRNwDfBnwH6tO9ML1fgHYpt32isA9JSm2ChT79wlwNXAbnYQxv9890wvX3O94D3R3vQs3QOwKjEhdDExv9x1/LSJuJBvu02/DmN2U3YiI5dptfh14spPTHgZWj4iVI2IBYB/AO6B6KKX0REppmZRSQ0qpgewvyI06BrGIWCIiBhdeLwUMB7xZooeKvd7AncCOheu+BFlL2p19XG5ViIjV223uBjzdyTl+v3tRMdccv+O9JiJ2Ar4L7JZSauvinIUjYtFZr8mud2f/tvYbhrHunVm4/fZxYFvgOICIWD4ibgcoDFY8muw/3qeACSmlf+RVcDWKiMaIuKiwuRYwOSIeA+4GzvDO1d7V/nqnlN4i66J/uPD4UWGfeu6MiHiy8PfJjsC3wO93iXV7zf2O96rzgEX5/+3dTahVZRSH8eevEFYSYZQ0iJSizAg/Iqis0LKQKAgzJIpqYEGTjJAgaBBIIweBA1H6gCIQErkEKpWYGIZm5bclNSihgRUOBIWwbDXYb3ivKNdEPafr8xtt1j577bMPh30W7/vus7qpx51JlsPQ30xgPLC5fce3AWur6pPevN3+4F9bSJIk9ZAjY5IkST1kMSZJktRDFmOSJEk9ZDEmSZLUQxZjkiRJPWQxJqmnkhw5R3kmtUfpdyS54VzklKQLwWJM0kjxGPBxVU0b3DooHe91kvqWNyhJfaEVTUvaH3TuSTK/xUclWZZkX5I1SdYlmXfSsQ8DLwMLkmxMMiHJ90mWAduB65I8lGRLku1JViUZ246dk2R/ks1JliZZ0+JvJFk06Bx7k0xo208n2dZG4lYkGd3iR1pT6l2tyff4Fh+fZKDFdyW5O8niJAsH5X8zyUvn7xOW1K8sxiT1i7nAVGAKMBtY0tqRzQUmALcBC4C7Tj6wqtYBy4G3qmpWC98MfFBV04CjwOvA7KqaDnwDvJJkDPA28ChwL6dufjxEkluA+XSNjqcCx4Gn2u7Lga1VNYWuz97zLb4U2NTi04F9wLvAsy3nKLo2aqdrqixpBLNRuKR+cQ+wsqqOA78m2QTc0eKrqupv4GCSjWeY70BVbW3bdwKTgS9bU/RLgC3AJOCnqvoRIMmHwAvD5H0AuB34uuW6FPit7TsGrGnb3wIPtu37gWcA2vUdBg4nOZRkGl17mB1VdegMr03SCGIxJqlf5D/Gh3P0pBzrq+rJIYmTqcDpesL9xdDZgzGDcr1fVa+d4pg/60SPueMMf499B3iObkTuvWFeK2mEcppSUr/4ApifZHSSq4H76JoIbwYeb2vHxgMzzyL3VmBGkhsBklyW5CZgPzBx0NOXg4u1n+mmFEkyHZjY4huAeUmuafvGJbl+mPNvAF5srx+d5IoWHwDm0I0AfnoW1yVpBLAYk9QvBoDdwC7gc+DVqjoIrAZ+AfYCK4Cv6Kb5zlhV/U43ArUyyW664mxSVf1BNy25Nslm4MCgw1YD45LspCukfmi5vqNbf/ZZy7UeuHaYt7AQmJVkD9305a0t1zFgI/BRm76UdBHKiRF1SepPScZW1ZEkV9GNls1ohdq5Ps9MYFFVPXKuc5/mfKPonvZ84t91a5IuPq4Zk/R/sCbJlXQL7xefj3pr4+QAAABCSURBVELsQksymW6x/4CFmHRxc2RMkiSph1wzJkmS1EMWY5IkST1kMSZJktRDFmOSJEk9ZDEmSZLUQxZjkiRJPfQPamoWm0Xag6oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#周期图\n",
    "\n",
    "per_list=[None for i in range(500)]\n",
    "\n",
    "for a in range(500):\n",
    "    nf = N/2\n",
    "    df = 1/(dt*N)\n",
    "    F_a = np.arange(1,nf+1)\n",
    "    F = [i*df for i in F_a]\n",
    "    F1 = F[0:int(nf)]\n",
    "    mean_x = np.mean(counts_list[a])\n",
    "    x_new  = [i-mean_x for i in counts_list[a]]\n",
    "    dft   = fft(counts_list[a])\n",
    "    dft_1 = dft[1:int(nf)+1]\n",
    "    per = (abs(dft_1)**2)\n",
    "\n",
    "    per_list[a] = per\n",
    "\n",
    "\n",
    "    \n",
    "# 500条周期图取 mean,std\n",
    "per_everypoint_list=[]\n",
    "for i in range(int(nf)):\n",
    "    per_everypoint_list.append([])\n",
    "ADM=[]\n",
    "ADM_std=[]\n",
    "for m in range(int(nf)):\n",
    "    for n in range(500):\n",
    "        per_everypoint_list[m].append(per_list[n][m])\n",
    "        \n",
    "for m in range(int(nf)):\n",
    "    ADM.append(np.mean(per_everypoint_list[m]))\n",
    "    ADM_std.append(np.std(per_everypoint_list[m]))    \n",
    "\n",
    "    \n",
    "\n",
    "\n",
    "POW1=POW[0:int(nf)]\n",
    "P_TIMES_F = np.multiply(np.array(F1),np.array(POW1))\n",
    "F1_log = [math.log(i,10) for i in F1]\n",
    "PTF_log = [math.log(i,10) for i in P_TIMES_F]\n",
    "\n",
    "\n",
    "\n",
    "# ADM取bin\n",
    "ADM_binned=databin_20(ADM)\n",
    "ADM_binned_std=databin_20_std(ADM)\n",
    "F1_binned=databin_20(F1)\n",
    "'''ADM_times_f_binned=np.multiply(np.array(F1_binned),np.array(ADM_binned))'''\n",
    "ADM_times_f_binned = np.array(F1_binned)+np.array(ADM_binned)\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(F1_log,PTF_log,color=\"b\",linewidth=1)\n",
    "plt.scatter(F1_binned,ADM_times_f_binned,color=\"r\",linewidth=1)    \n",
    "plt.xlabel(\"log frequency\")\n",
    "plt.ylabel(\"log power*frequency\")\n",
    "plt.title(\"power density spectrum\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAHwCAYAAADNfOnlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xm8nOP9//HXJ5slQpDYQhJrN/yCY1e7SLWSVGzJQVQISmmLqsau+VZp8aX6JW1qTSyxNai9aislUWKpWlqJohV7KpZErt8f9xwZx1kmOWfOPXPm9Xw85nFm7vueM+9Jinfv676uO1JKSJIkqTp1yTuAJEmSFp9lTpIkqYpZ5iRJkqqYZU6SJKmKWeYkSZKqmGVOkiSpilnmJHVKEXFpRPy0zJ/xTERsX87PkKTWWOYkaTGllL6WUvoTQEScGhFX5hypZBHxckTsnHcOSW1nmZNUVSLjv7vKLCK65Z1BUmn8F6KkVhXO4pwQEc9GxDsRcUlELFm0/5CIeDEi3o6IqRGxWmH7aRFxQeF594j4ICLOKrxeKiI+iojlC6+3iIg/R8S7EfFk8fBlRPwpIsZHxEPAXGCtJjJuFBGPR8SciLgGWLLR/m9FxBOF3//niNiw0fc7NiJmRMR7EXFNw/eLiD4RcUvhfW9HxAMNZbLh7FZEDAF+AuwTEf8t5N8rIqY3ynBMRNzUzJ/xgRHxj0L+f0ZEfdH2hyLigkK25yJip6L3LRcREyPi9Yh4NSJ+GhFdG/3d/K3we5+NiI0j4gqgP3BzIe+PImJgRKSIGBMRs4A/RsT2EfGvJv63sHPh+akRMSUiriz8/qciYr3C/1beiIhXImJwU99XUvuxzEkqVT2wK7A2sB5wIkBE7Aj8DNgbWBWYCVxdeM99wPaF55sC/wa2K7zeEvh7SumdiOgH3Ar8FFgBOBa4PiL6Fn3+/sBYoFfhMz4TET2Am4ArCu+fAowo2r8x8DvgUGBF4GJgakQsUfRr9gaGAGsCGwIHFrYfA/wL6AusTFbaPncfxJTS7cD/ANeklJZJKf0/YCqwZkR8pejQ/QoZPyciegLnA99IKfUCtgKeKDpkc+AfQB/gFOCGiFihsO8yYD6wDrARMBg4uPB79wJOBQ4AlgWGAm+llPYHZgG7F/KeVfRZ2wFfIfu7LsXuhe+0PPBX4A6y/7b0A04n+7OWVEaWOUml+lVK6ZWU0tvAeGBkYXs98LuU0uMppY+BE4AtI2Ig8DCwbkSsCGwLTAT6RcQyZKXhvsLv2A/4Q0rpDymlBSmlu4BpwG5Fn39pSumZlNL8lNK8Rtm2ALoD56WU5qWUrgMeK9p/CHBxSukvKaVPU0qXAR8X3tfg/JTSa4XvdzMwqLB9HllJHVD43Q+kEm5qXfizuKbw3YiIrwEDgVuaecsCYP2IWCql9HpK6ZmifW8UfbdrgL8D34yIlYFvAN9PKX2QUnoDOBfYt/C+g4GzUkqPpcyLKaXPFeEmnFr4XR+29h0LHkgp3ZFSmk9WovsCZxb+jq4GBkZE7xJ/l6TFYJmTVKpXip7PBFYrPF+NojNlKaX/Am8B/QqFYBpZcduWrLz9Gdiaz5e5AcBehaHMdyPiXWAbshLV1Oc3thrwaqOSVVxaBgDHNPr9axR9B8jOGjaYCyxTeH428CJwZ2EY9Mct5GjsMmBURATZmcVrCyXvc1JKHwD7AIcBr0fErRHx5aJDmvpuqxW+V/fCexq+18XASoXj1gBeWoS80PKfc1P+U/T8Q+DNlNKnRa9h4Z+lpDKwzEkq1RpFz/sDrxWev0ZWKoDPhgxXBF4tbLoP2JFsCPCxwutdgc2A+wvHvAJckVLqXfTomVI6s+gzWzob9jrZGb9olLHBK8D4Rr9/6ZTSVa196ZTSnJTSMSmltciGFH9YfM1aS/lSSo8AnwBfB0bRxBBr0bF3pJR2ISuwzwG/Kdrd1Hd7rfC9Pgb6FH2vZVNKXyv63ms395ElbP8AWLrhReFavL5feIekXFnmJJXqiIhYvXCt1k/IhhABJgPfiYhBhWvQ/gf4S0rp5cL++8iu2Xo2pfQJ8Cey4b9/ppRmF465Etg9InaNiK4RsWTh4vvVS8z2MNl1Y0dFRLeI2IOsLDb4DXBYRGwemZ4R8c2I6NXaL45s4sQ6hTL1PvBp4dHYf8iGFBv/e/Vy4FfA/JTSg818xsoRMbRQhD8G/tvoM1YqfLfuhevgvkI2LP06cCfwy4hYNiK6RMTaEdFwXeJvgWMjYpPC914nIhqK939oYiJJI88DSxb+rLqTXSe5RCvvkdTBLHOSSjWZrDj8o/D4KUBK6R7gJOB6sjNka7Pwmi3IhlWXYuFZuGeBj4pek1J6BRhGVhJnk51ROo4S/x1VKIl7kE1aeIdsyPKGov3TyK6b+1Vh/4ssnODQmnWBu8kK1sPArxvWlmtkSuHnWxHxeNH2K4D1aeGsHNn3PIbsbNvbZEPQ3y3a/5dCjjfJrlfcM6X0VmHfAUAPsj/Xd4DrKAxPp5SmFI6fDMwhmyTSMHHiZ8CJheHZY5sKlVJ6r5Djt2RnWj8gmwwiqYJECdfxSqpxEfEycHBK6e68s1SbiFiKbALDximlFxbj/QeS/dlv097ZJHUOnpmTpPI6HHhscYqcJJXCFb4lqUwKZzQDGJ5zFEmdmMOskiRJVcxhVkmSpCpmmZMkSapiNXXNXJ8+fdLAgQPzjiFJktSq6dOnv5lSanWh7poqcwMHDmTatGl5x5AkSWpVRLR2L2XAYVZJkqSqZpmTJEmqYpY5SZKkKmaZkyRJqmKWOUmSpCpmmZMkSapiljlJkqQqZpmTJEmqYpY5SZKkKmaZkyRJqmKWOUmSpCpmmZMkSapiljlJkqQqZpmTJEmqYpY5SZKkKmaZkyRJqmKWOUmSpCpmmZMkSapiljmplkyaBAMHQpcu2c9Jk/JOJElqo255B5DUQSZNgrFjYe7c7PXMmdlrgPr6/HJJktrEM3NSrRg3bmGRazB3brZdklS1LHNSrZg1a9G2S5KqQi5lLiLOjojnImJGRNwYEb2bOGaNiLg3Iv4WEc9ExNFF+06NiFcj4onCY7eO/QZSFerff9G2S5KqQl5n5u4C1k8pbQg8D5zQxDHzgWNSSl8BtgCOiIivFu0/N6U0qPD4Q/kjS1Vu/HhYeunPb1t66Wy7JKlq5VLmUkp3ppTmF14+AqzexDGvp5QeLzyfA/wN6NdxKaVOpr4eJkyAAQMgIvs5YYKTHySpylXCbNaDgGtaOiAiBgIbAX8p2nxkRBwATCM7g/dOuQJKnUZ9veVNkjqZsp2Zi4i7I+LpJh7Dio4ZRzac2uxiVxGxDHA98P2U0vuFzf8HrA0MAl4HftnC+8dGxLSImDZ79ux2+GaSJEmVI1JK+XxwxGjgMGCnlNLcZo7pDtwC3JFSOqeZYwYCt6SU1m/tM+vq6tK0adMWO7MkSVJHiYjpKaW61o7LazbrEOB4YGgLRS6AicDfGhe5iFi16OW3gafLlVWSJKmS5TWb9VdAL+CuwtIiFwFExGoR0TAzdWtgf2DHJpYgOSsinoqIGcAOwA86+gtIkiRVglwmQKSU1mlm+2vAboXnDwLRzHH7ly+dJElS9fAOEJIkSVXMMidJklTFLHOSJElVzDInSZJUxSxzkiRJVawSbufVacybB2++mT2PyB5ted7W9zd+LkmSOh/LXDv6+99hl10gpewBi/+8re9veN5YWwtjly7QtWt+j27dYMklFz6WWKL51y3tW3LJ7HdZdCVJ1c4y147WXx9efz3vFE1rr8K4YAF8+ml+j3nz4OOPs8f778NHH2WPjz9e+Lzx6+b2ffrp58vdUktB796w3HLZz1KfL7tsVjQlScqDZa5GONz6RfPnLyyGH30EH3wA772XPd59N3s0PH/ppaa3v/suzJkDPXsuLHl9+kC/frD66tnP4scqq2RnBCVJai/+Z0U1q1u37NGzZ9t+z4IFWaFrKHmzZ8Orr2aP556De+5Z+PrNN6Fv3y+WvIby96UvwWqrWbwlSaWzzElt1KVLdkZuueVaP3b+fPj3vxeWu3/9K/v5zDPZ8+eey84SfvWr2eNrX8seX/1qVvgseZKkxixzUgfq1i07A7f66s0f89Zb8OyzWcF75hm45Zbs54cffrHgbbxxdqZPklS7LHNShVlxRfj617NHscYl7+ab4fHHs2HZbbeF7bbLfvbrl09uSVI+LHNSlWiq5H36KTz5JNx/P0yZAt/7Xjbc21Dstt0W1lzT4VlJ6swiNbUYWSdVV1eXpk2blncMqWwWLIC//S0rd/fdl/3s2jUrdUOHwje/Ccssk3dKSVIpImJ6Sqmu1eMsc1LnlVK2rMq998INN8Cf/wyDB8Pee8Nuu7V9Jq8kqXxKLXPem1XqxCJgnXXgkEPgttvgn/+Eb3wDJk7MrrXbe2+47jqYOzfvpJKkxWWZk2rICivAQQfB7bdnZ+wGD4YJE7Jit+++2czZBQvyTilJWhSWOalG9ekDBx8Md94JL74IO+4Ip50GX/4yXHhhdkcMSVLls8xJok8fGDsWHn0Ufvc7+OMfYcAAOP74bDFjSVLlssxJ+kwEbLMNXH99Vuw+/hg23BBGjYLHHss7nSSpKZY5SU1aay0477xs0kRdHey5Z7bG3QMP5J1MklTMMiepRcstBz/8YTZh4tBDYb/9smL30kt5J5MkgWVOUom6dcuK3HPPZfeE3XxzOPZYePfdvJNJUm2zzElaJEstBT/5CTz9NMyZA1/6EvzqVzBvXt7JJKk2WeYkLZZVVoGLL4a774apU2GDDbKFiSVJHcsyJ6lNNtgA7rgDzjkHvvvd7G4Tc+bknUqSaodlTlKbRWT3en3ySfj0Uxg0CB56KO9UklQbLHOS2s2yy2aLDp9zTjbj9YQT4JNP8k4lSZ2bZU5Suxs2LDtL9+yzsNlm2WQJSVJ5WOYklcVKK8FNN8H3vgc77ADnnw8p5Z1Kkjofy5yksomAMWPgL3+BiROzyREOu0pS+7LMSSq7tdbKJkS89RbsvDPMnp13IknqPCxzkjrEMsvA9ddn93fdbDN46qm8E0lS52CZk9RhunSB8eOzx047ZYsNS5LaplveASTVnlGjYO21YY894G9/gx/9KLu+TpK06DwzJykXm2+eTYyYPBmOO86ZrpK0uCxzknKz+upw773wpz/B0Udb6CRpcVjmJOVqhRXgnnvgscfgsMNgwYK8E0lSdbHMScrdcsvBnXdm188ddFB2f1dJUmksc5IqQq9ecNtt8MorsP/+MH9+3okkqTpY5iRVjJ494ZZb4J13YN99LXSSVArLnKSKstRS2T1d58yBww93UoQktcYyJ6niLLEEXHcd/PWvcNppeaeRpMrmosGSKlKvXnDrrbD11rDaajB2bN6JJKkyWeYkVayVV4bbb8/u57ryyjBsWN6JJKnyWOYkVbR11snu4brbbtC3L2y1Vd6JJKmyeM2cpIq36aZwxRUL7+UqSVoolzIXEWdHxHMRMSMiboyI3s0c93JEPBURT0TEtKLtK0TEXRHxQuHn8h2XXlIehgyBn/0sG2p9992800hS5cjrzNxdwPoppQ2B54ETWjh2h5TSoJRSXdG2HwP3pJTWBe4pvJbUyX3nOzB4cLaosLf9kqRMLmUupXRnSqlhOdBHgNUX8VcMAy4rPL8MGN5e2SRVtnPOyRYVPuOMvJNIUmWohGvmDgJua2ZfAu6MiOkRUbwwwcoppdcBCj9Xau6XR8TYiJgWEdNmz57dbqEl5aNHj2wNut/8JrtbhCTVurKVuYi4OyKebuIxrOiYccB8YFIzv2brlNLGwDeAIyJi20XNkVKakFKqSynV9e3bd7G+i6TKssoqMGUKHHQQvPBC3mkkKV9lW5okpbRzS/sjYjTwLWCnlJq+YU9K6bXCzzci4kZgM+B+4D8RsWpK6fWIWBV4o33TS6p0W24Jp58Ow4fDX/4CyyyTdyJJykdes1mHAMcDQ1NKc5s5pmdE9Gp4DgwGni7sngqMLjwfDfy+vIklVaJDD81K3SGHeA9XSbUrr2vmfgX0Au4qLDtyEUBErBYRfygcszLwYEQ8CTwK3JpSur2w70xgl4h4Adil8FpSjYmACy6Ap56Cyy/PO40k5SOaGeHslOrq6tK0adNaP1BSVZkxA3baCR55BNZeO+80ktQ+ImJ6o6XZmlQJs1klqU023BBOPBHq62HevLzTSFLHssxJ6hS+9z3o3dv15yTVHsucpE6hSxe49NJs/bkHHsg7jSR1HMucpE5jlVWyMrf//t6/VVLtsMxJ6lS+9S345jfh6KPzTiJJHcMyJ6nT+fnP4f774fbbWz9WkqqdZU5Sp7PMMjBhQrao8Jw5eaeRpPKyzEnqlHbZBXbeGU44Ie8kklReljlJndYvfgE33ujsVkmdm2VOUqe1/PJw4YUwZgx8+GHeaSSpPCxzkjq14cNh0CA47bS8k0hSeVjmJHV6F1wAl1wCTz6ZdxJJan+WOUmd3sorw+mnw5FHQkp5p5Gk9mWZk1QTDj44u27uyivzTiJJ7csyJ6kmdO0Kv/41HH88vPde3mkkqf1Y5iTVjM02y271dcopeSeRpPZjmZNUU372M5g8GWbMyDuJJLUPy5ykmtKnTzYZ4ogjnAwhqXOwzEmqOYcckk2GmDw57ySS1HaWOUk1p2tXOO+87L6t3hlCUrWzzEmqSdtsA5tumpU6SapmljlJNevnP4df/hLeeCPvJJK0+CxzkmrWOuvA/vu7VImk6maZk1TTTjoJrr8enn027ySStHgsc5Jq2gorZBMhfvSjvJNI0uKxzEmqeUccAc89B/fck3cSSVp0ljlJNa9HD/if/8nu2+pCwpKqjWVOkoA994QFC+DGG/NOIkmLxjInSUCXLjB+PJx4Inz6ad5pJKl0ljlJKhgyBFZcESZNyjuJJJXOMidJBRHZtXOnnAKffJJ3GkkqjWVOkop8/evw5S/Db36TdxJJKo1lTpIa+elPs+vnPvgg7ySS1DrLnCQ1sskmsPXWcOGFeSeRpNZZ5iSpCaecAuec49k5SZXPMidJTVh//ez6uYsuyjuJJLXMMidJzTjpJPjFL2Du3LyTSFLzLHOS1IwNN4Qtt4SLL847iSQ1zzInSS04+WQ4+2z48MO8k0hS0yxzktSCQYNgs81gwoS8k0hS0yxzktSKk0+Gs86Cjz7KO4kkfZFlTpJasfHGsNFGcOmleSeRpC+yzElSCU44ITs7N39+3kkk6fMsc5JUgq23htVXh2uvzTuJJH2eZU6SSnTCCXDmmZBS3kkkaSHLnCSVaMgQ6NoVbr017ySStJBlTpJKFAE//nF2dk6SKkUuZS4izo6I5yJiRkTcGBG9mzjmSxHxRNHj/Yj4fmHfqRHxatG+3Tr+W0iqRXvuCf/+NzzwQN5JJCmT15m5u4D1U0obAs8DJzQ+IKX095TSoJTSIGATYC5wY9Eh5zbsTyn9oUNSS6p5XbvCj37k2TlJlSOXMpdSujOl1DDB/xFg9VbeshPwUkppZnmTSVLrDjgApk+HZ57JO4kkVcY1cwcBt7VyzL7AVY22HVkYpv1dRCxfnmiS9EVLLglHHgnnnJN3EkmCSGWaYx8RdwOrNLFrXErp94VjxgF1wB6pmSAR0QN4DfhaSuk/hW0rA28CCTgDWDWldFAz7x8LjAXo37//JjNnenJPUtu99Rasu252dm7VVfNOI6kziojpKaW6Vo8rV5lr9YMjRgOHATullOa2cNww4IiU0uBm9g8Ebkkprd/aZ9bV1aVp06YtXmBJauTII2G55WD8+LyTSOqMSi1zec1mHQIcDwxtqcgVjKTREGtEFP//4G8DT7dvQklq3Q9+ABMmwAcf5J1EUi3L65q5XwG9gLsKS4tcBBARq0XEZzNTI2JpYBfghkbvPysinoqIGcAOwA86KLckfWbttWG77eCSS/JOIqmW5TbMmgeHWSW1t0cegVGj4PnnoVu3vNNI6kwqephVkjqLLbaAVVaBm2/OO4mkWmWZk6Q2OuoouOCCvFNIqlWWOUlqoxEj4O9/h6eeyjuJpFpkmZOkNureHQ47zLNzkvJhmZOkdjB2LEyZki0mLEkdyTInSe1g5ZVh6FCYODHvJJJqjWVOktrJUUfBr38N8+fnnURSLbHMSVI72WQTWG01lymR1LEsc5LUjo46Cs4/P+8UkmqJZU6S2tGIEdndIGbMyDuJpFphmZOkdtS9Oxx+uMuUSOo4ljlJamdjx8J117lMiaSOYZmTpHa20kouUyKp41jmJKkMjjoKLrzQZUoklZ9lTpLKYJNNoF8/+P3v804iqbOzzElSmRx5JPzf/+WdQlJnZ5mTpDIZMSJbouSFF/JOIqkzs8xJUpkssQQceCBMmJB3EkmdmWVOkspo7Fi49FL46KO8k0jqrCxzklRG66wDgwbBDTfknURSZ9VqmYuIFToiiCR1VoceChddlHcKSZ1VKWfm/hIRUyJit4iIsieSpE5m2LBsEsSzz+adRFJnVEqZWw+YAOwPvBgR/xMR65U3liR1Ht27w5gxcPHFeSeR1Bm1WuZS5q6U0kjgYGA08GhE3BcRW5Y9oSR1AoccAldeCXPn5p1EUmdTyjVzK0bE0RExDTgW+B7QBzgGmFzmfJLUKQwYAFtsAddem3cSSZ1NKcOsDwPLAsNTSt9MKd2QUpqfUpoGeEmvJJXo0EMdapXU/kopc19KKZ2RUvpX4x0ppZ+XIZMkdUq77Qb/+hc8+WTeSSR1JqWUuTsjonfDi4hYPiLuKGMmSeqUunWDgw/27Jyk9lVKmeubUnq34UVK6R1gpfJFkqTO6+CD4eqrYc6cvJNI6ixKKXOfRkT/hhcRMQBI5YskSZ1Xv36w3XZw1VV5J5HUWZRS5sYBD0bEFRFxBXA/cEJ5Y0lS53XwwTBxYt4pJHUWpawzdzuwMXANcC2wSUrJa+YkaTHtums2EeLpp/NOIqkzKOXMHMASwNvAe8BXI2Lb8kWSpM6tWzc48EDPzklqH6UsGvxz4CGy4dbjCo9jy5xLkjq1gw7K7gjx8cfApEkwcCB06ZL9nDQp53SSqkm3Eo4ZTrbW3MflDiNJtWLttWGDDWDqcQ+w18SxC+/zNXMmjB2bPa+vzy+gpKpRyjDrP4Du5Q4iSbVmzBiY+JsFX7xh69y5MG5cPqEkVZ1SzszNBZ6IiHuAz87OpZSOKlsqSaoBe+wBR+23PrNYg/688vmds2blE0pS1SnlzNxU4Azgz8D0oockqQ2WWgr2XeZWLuXAL+7s3/+L2ySpCa2emUspXRYRSwH9U0p/74BMklQzxhzfhxEnb8uJ6ad0aViPfemlYfz4fINJqhqlzGbdHXgCuL3welBETC13MEmqBRufuBu9+y/LH1caCREwYABMmODkB0klK2WY9VRgM+BdgJTSE8CaZcwkSTVlzLErMHHHSbBgAbz8skVO0iIppczNTym912ib92aVpHZSXw+33QZvvZV3EknVqJQy93REjAK6RsS6EXEB2WQISVI7WH552G031wqWtHhKKXPfA75GtizJVcD7wPfLGUqSas2YMdntvZLjHpIWUSmzWeeS3crLFSwlqUx22AHefx+eeAI22ijvNJKqSatlLiLupYlr5FJKO5YlkSTVoC5dYP/94bLLLHOSFk0pd4A4tuj5ksAIYH554khS7TrgANh6azj7bOjuTRQllaiUYdbGd3t4KCLuK1MeSapZ66wD666bzWwdOjTvNJKqRSmLBq9Q9OgTEbsCq7T1gyPijIiYERFPRMSdEbFaM8eNjogXCo/RRds3iYinIuLFiDg/IqKtmSQpbwcemA21SlKpIrUydSoi/kl2zVyQDa/+Ezg9pfRgmz44YtmU0vuF50cBX00pHdbomBWAaUBdIcN0YJOU0jsR8ShwNPAI8Afg/JTSbS19Zl1dXZo2bVpbYktSWb33XnYTiJdeghVXzDuNpDxFxPSUUl1rx7V6Zi6ltGZKaa3Cz3VTSoPbWuQKv/f9opc9aXoh4l2Bu1JKb6eU3gHuAoZExKrAsimlh1PWRi8Hhrc1kyTlbbnlsjXnrroq7ySSqkUps1n3aGl/SumGxf3wiBgPHAC8B+zQxCH9gFeKXv+rsK1f4Xnj7U19xlhgLED//v0XN6okdZjRo+HEE+HII/NOIqkalLJo8BhgIlBfePwW2A/YHfhWS2+MiLsj4ukmHsMAUkrjUkprAJOApv611dR1cKmF7V/cmNKElFJdSqmub9++LcWVpIqw887w2mvw7LN5J5FUDUpZmiSRXc/2OkBhiPPClNJ3Wn1jSjuXmGMycCtwSqPt/wK2L3q9OvCnwvbVG21/rcTPkqSK1rUr7LdfNhHi5z/PO42kSlfKmbmBDUWu4D/Aem394IhYt+jlUOC5Jg67AxgcEctHxPLAYOCOQp45EbFFYRbrAcDv25pJkirF6NFw5ZXw6ad5J5FU6Uopc3+KiDsi4sDC0iC3Ave2w2efWRhynUFW0o4GiIi6iPgtQErpbeAM4LHC4/TCNoDDyYZ8XwReAlqcySpJ1eSrX4V+/eDuu/NOIqnStbo0CUBEfBvYtvDy/pTSjWVNVSYuTSKpmlx4ITz0EEyenHcSSXkodWmSUq6ZA3gcmJNSujsilo6IXimlOW2LKElqyb77wrhx2dpzyy2XdxpJlaqUO0AcAlwHXFzY1A+4qZyhJEnZosE77QTXXpt3EkmVrJRr5o4AtgbeB0gpvQCsVM5QkqTM6NHe3ktSy0opcx+nlD5peBER3WhmTTdJUvv6xjfghRfgxRfzTiKpUpVS5u6LiJ8AS0XELsAU4ObyxpIkAXTvDiNHwuWX551EUqUqpcz9GJgNPAUcSnZT+xPLGUqStNDo0VmZW7Ag7ySSKlGLs1kjoitwWUppP+A3HRNJklRs0KBsNuv998P22+edRlKlafHMXErpU6BvRPTooDySpEYinAghqXmlDLO+DDwUESdFxA8bHmXOJUkqsu++cNNN8OGHeSeRVGlKKXOvAbcUju1V9JAkdZDVVoNNN4WbnX4mqZFmr5mLiCtSSvsD76aU/rcDM0mSmrDffjBpEuy9d95JJFWSls7MbRIRA4CDImL5iFih+NFRASVJmW9oD+OrAAAgAElEQVR/G+67D958M+8kkipJS2XuIuB24Mtk92adXvTwbvWS1MF69coWEZ4yJe8kkipJs2UupXR+SukrwO9SSms2eqzVgRklSQUNQ62S1KDFCRARMSqldHhE7NtRgSRJzRs8GJ5/Hv7xj7yTSKoUrc1m7RcRewOrd0QYSVLLunfPJkBMnpx3EkmVotkyFxGnACsAk4EVIuLkDkslSWpWw1BrSnknkVQJWrpm7jTgbWA/4O2U0ukdlkqS1KzNN4d58+Dxx/NOIqkStDbM+lpK6Wrg1Y4II0lqXQTU18OVV+adRFIlaK3MXRURT6eUruqQNJKkktTXw9VXw/z5eSeRlLcWy1xKaQHwZET076A8kqQSrLcerLEG/PGPeSeRlLdmb+dVZFXgmYh4FPigYWNKaWjZUkmSWrXfftlQ6+DBeSeRlKdSytxpZU8hSVpk++wDJ58MH3wAPXvmnUZSXlq7Zo6U0n3Ay0D3wvPHyG7vJUnK0corw5ZbwtSpeSeRlKdWy1xEHAJcB1xc2NQPuKmcoSRJpWkYapVUu1otc8ARwNbA+wAppReAlcoZSpJUmuHD4aGH4I038k4iKS+llLmPU0qfNLyIiG6A645LUgXo2RO+9S249tq8k0jKSyll7r6I+AmwVETsAkwBbi5vLElSqRxqlWpbKWXux8Bs4CngUOAPwInlDCVJKt3OO8PLL8MLL+SdRFIeSilz2wOTUkp7pZT2TCn9JiVv7yxJlaJbN9h3X5g8Oe8kkvJQSpk7EHgiIh6OiLMiYveIWL7MuSRJi6DhXq3+X22p9pSyztwBKaX1gBHAv4ALyYZdJUkVoq4OunSBRx/NO4mkjtbqHSAiYj/g68AGwJvAr4AHypxLkrQIIhZOhNh887zTSOpIpdzO6zzgJeAi4N6U0stlTSRJWiyjRsFWW8G552bX0UmqDaUMs/YBDgKWBMZHxKMRcUXZk0mSFsnaa8Oaa8Ldd+edRFJHKuV2XssC/YEBwEBgOWBBeWNJkhZHfb2zWqVaU8ps1geB3YEZwD4ppS+llEaXN5YkaXHsvTfcfDPMnZt3EkkdpdWrKlJKGwJERC+8jZckVbSVV4bNNssK3T775J1GUkcoZZh1/Yj4K/A08GxETI+I9csfTZK0OBxqlWpLKcOsE4AfppQGpJT6A8cUtkmSKtDw4XDfffD223knkdQRSilzPVNK9za8SCn9CehZtkSSpDZZdlnYdVe47rq8k0jqCKWUuX9ExEkRMbDwOBH4Z7mDSZIW36hRDrVKtaKUMncQ0Be4Abix8Pw75QwlSWqbIUPgqafglVfyTiKp3EpZNPidlNJRwA7Atimlo1NK75Q/miRpcS2xBIwYAVdfnXcSSeVWymzWTSPiKeBJ4KmIeDIiNil/NElSW4waBZMm5Z1CUrmVMsw6EfhuSmlgSmkgcARwSVlTSZLabNtt4c034Zln8k4iqZxKKXNzUkoPNLxIKT0IzClfJElSe+jSBUaOdCKE1NmVUuYejYiLI2L7iNguIn4N/CkiNo6IjRfnQyPijIiYERFPRMSdEbFaE8cMioiHI+KZwrH7FO27NCL+WXj/ExExaHFySFJn1zCrNXn/HqnTitTKP+ERcW8Lu1NKacdF/tCIZVNK7xeeHwV8NaV0WKNj1iv8/hcKZW868JWU0rsRcSlwS0ppkVZRqqurS9OmTVvUuJJUtVKCr30NJk6ELbfMO42kRRER01NKda0dV8q9WXdon0if+53vF73sSRP3fE0pPV/0/LWIeINsWZR32zuPJHVWEQsnQljmpM6plGHWsoiI8RHxClAPnNzKsZsBPYCXijaPLwy/nhsRS5QxqiRVtVGj4NprYd68vJNIKoeylbmIuDsinm7iMQwgpTQupbQGMAk4soXfsypwBfCdlNKCwuYTgC8DmwIrAMe38P6xETEtIqbNnj27nb6dJFWPtdaCtdeGu+/OO4mkcmj1mrmyB4gYANyaUlq/iX3LAn8CfpZSmtLM+7cHjk0pfau1z/KaOUm16oIL4NFH4Yor8k4iqVTtds1cROzRxOb3gKdSSm8sZrh1U0ovFF4OBZ5r4pgeZLcPu7xxkYuIVVNKr0dEAMOBpxcnhyTVir33hpNOgrlzYeml804jqT21WuaAMcCWQMOs1u2BR4D1IuL0lNLi/P+8MyPiS8ACYCZwGEBE1AGHpZQOBvYGtgVWjIgDC+87MKX0BDApIvoCATzR8H5JUtNWXhm22AKmToV99807jaT2VEqZW0C2JMh/ACJiZeD/gM2B+8muZ1skKaURzWyfBhxceH4lcGUzxy3yciiSVOvq67M15yxzUudSygSIgQ1FruANYL2U0tuAc6MkqUoMHw733QdvvZV3EkntqZQy90BE3BIRoyNiNDAVuD8ieuKab5JUNXr1giFD4LpFWm5dUqUrpcwdAVwCDAI2Ai4DjkgpfVCOBYUlSeXTMNQqqfNotcylbO2SB4E/AncD96e81zORJC2WIUPgmWdg1qy8k0hqL62WuYjYG3gU2JNshulfImLPcgeTJLW/Hj1gxAi4+uq8k0hqL6UMs44DNk0pjU4pHQBsBpxU3liSpHJpuFerpM6hlDLXpdHiwG+V+D5JUgX6+tfh7bfhaZdblzqFUkrZ7RFxR0QcWFi891bgD+WNJUkqly5dYORIuOqqvJNIag+lTIA4DpgAbAj8P2BCSqnZG9tLkirfqFHZrFans0nVr5Q7QJBSuh64vsxZJEkd5P/9P1hqKXj4Ydhqq7zTSGqLZs/MRcSciHi/iceciHi/I0NKktpXRLbmnBMhpOrXbJlLKfVKKS3bxKNXSmnZjgwpSWp/I0fClCkwzxszSlXNWamSVKPWWgvWWQfuvjvvJJLawjInSTXMNeek6meZk6QatvfecMst8MEHeSeRtLgsc5JUw1ZaCbbcEqZOzTuJpMVlmZOkGtew5pyk6mSZk6QaN3w4PPAAvPVW3kkkLQ7LnCTVuF69YMiQbJkSSdXHMidJcqhVqmKWOUkSQ4bAs8/CzJl5J5G0qCxzkiR69IARI+Dqq/NOImlRWeYkSUB2r1aHWqXqY5mTJAGwzTbw9tvw9NN5J5G0KCxzkiQAunSBkSM9OydVG8ucJOkzDUOtCxbknURSqSxzkqTPbLgh9OwJDz+cdxJJpbLMSZI+E+FECKnaWOYkSZ8zcmR2N4h58/JOIqkUljlJ0uesuSassw7cdVfeSSSVwjInSfqC+nqYNCnvFJJKYZmTJH3BXnvBrbfCBx/knURSayxzkqQvWGkl2GormDo17ySSWmOZkyQ1adQoh1olIPsHYeDAbGXtgQMr7h8My5wkqUnDhsEDD8Cbb+adRMrRpEkwdixvzvwv49MJpJkzYezYiip0ljlJUpN69YLddsuWKZFq1rhxzJv7CXsxhTn0IgDmzoVx4/JO9hnLnCSpWaNGuYCwalTD0OrMmXyf8+jJB4ynqMDNmpVbtMYsc5KkZu26K/ztbzBzZt5JpA5UGFpl5kwuZiz3sgOTGUVXim5a3L9/fvkascxJkprVowfsuSdcdVXeSaQONG4czJ3L/XydkzmdqQxlWeYs3L/00jB+fH75GrHMSZJa5FCrakbR0OrLDGAfruFK9mMdXlp4zIABMGFCtrJ2heiWdwBJUmXbZht491146inYYIO800hl0jC0Oncu/6Unw/g9P+ZMduHuhccMGAAvv5xbxOZ4Zk6S1KIuXWDkSJg87pmKXmtLapPC0OoCgv25gk15jKM4f+H+ChtaLWaZkyS1qr73rUy+pRcLZs6ClLIZERW21pbUJoXZqadyKm/Sh1/z3WwZEqjIodViljlJUqs2uOgIeqX3+TNbLdxYYWttSYul4Tq5lLiWvbicA7ieEfRgXra/YWi1QoscWOYkSSWIV2YxislMZtTnd1TQWlvSIitaguRxNuIILuQmhrMSs7P9FTy0WswyJ0lqXf/+jOQqprAX84rnzlXQWlvSIitcJ/dvVmY4N3ERhzGIJ7N9FT60WswyJ0lq3fjxrLn0G6zH89zJ4GxblZy1kJo1axYf04M9uIExTGQEN2TbIyp+aLWYZU6S1Lr6epgwgVEr3M4k6qvqrIX0BYXr5FJKHMZF9ONVTuKMhfur7IxzbuvMRcQZwDBgAfAGcGBK6bUmjvsUeKrwclZKaWhh+5rA1cAKwOPA/imlTzoiuyTVpPp69h4MP1kH/vv0KJZZJu9A0mIoWk/uXH7AEwziQbahCynbX4VnnPM8M3d2SmnDlNIg4Bbg5GaO+zClNKjwGFq0/efAuSmldYF3gDFlzitJNa9vX9h6a5g6Ne8k0mIqXCd3G0P4Bcfye4bRk7nZvio945xbmUspvV/0sic0VOLWRUQAOwLXFTZdBgxvv3SSpObU17u8nKrYrFk8x5cYzWVMYS/680q2vcqukyuW6zVzETE+Il4B6mn+zNySETEtIh6JiIbCtiLwbkppfuH1v4B+zXzG2ML7p82ePbtd80tSLRo2DB58EPxXqqpGw1pyXbrwTqzAUKZyJj9ma/688Jgqu06uWFnLXETcHRFPN/EYBpBSGpdSWgOYBBzZzK/pn1KqA0YB50XE2rBwUeYiTZ7ZSylNSCnVpZTq+vbt2w7fSpJq2zLLwG67wZQpeSeRSlC0ltz81IV9Fkzmm9zKQVyy8JgqvE6uWFnLXEpp55TS+k08ft/o0MnAiGZ+x2uFn/8A/gRsBLwJ9I6IhgkcqwNfmDwhSSoPh1pVNQrXyAEcx9kAnM1x0LVrNrRapdfJFcttmDUi1i16ORR4roljlo+IJQrP+wBbA8+mlBJwL7Bn4dDRQOOCKEkqk8GD4YUX4KWX8k4itaJwl5Lf8R3+wG5cwz5041NYsCB7VOl1csXyvGbuzMKQ6wxgMHA0QETURcRvC8d8BZgWEU+SlbczU0rPFvYdD/wwIl4ku4ZuYsfGl6Ta1aMHjBwJl1+edxKpFf378xBb8WPOZCpDWZ53P9veWUR2kqs21NXVpWnTpuUdQ5I6hccfhxEjsrNzXVyCXpVm0iQYN45ZMxewBY/wOw5iCHdk+5ZeuiqGViNiemHeQIv8x0+StFg22iibDPHAA3knkRopTHr4YOZshjKVY/kFQ+LObF8nuEausdzuACFJqm4RcOCBcOmlsN12eaeRiowbx4K5HzKaa9mIv/IDzs3WvBgwILtGrpPxzJwkabHV18NNN8F//5t3EqnIrFmcwUm8xmpcxGEL1zMrTIbobCxzkqTFtsoq2e29brgh7ySqeUULA18fI5jIGG5gD5ag6LbtnWjSQzHLnCSpTUaPhssuyzuFalrRwsBPpA05bMGvuYnhrMJ/Fh5T5QsDt8QyJ0lqk913hyefhJkz806imlVYGPgN+jKcm7iQI9iYv3aqhYFbYpmTJLXJkkvC3nvDFVfknUQ1a9YsPqE7e3ADB3A5e1O411wnWhi4JZY5SVKbNQy11tDSpaogaY3+fJdfsxJvcCqnLtzRSa+Ra8wyJ0lqs802g27d4M9/zjuJakbRhIfz3xzFY7EZl3MAXSj8P4pOfI1cY5Y5SVKbRTgRQh2oaMLDnWlnzpz7PX7fdQTLrLhkTVwj15iLBkuS2sX++8MGG8D//i8stVTeadSpFSY8PM+67M8VTGEvBs5/EZYZAG++mXe6DueZOUlSu+jXDzbdNFtEWCqrWbN4l+UYylR+yolsywOfba9FljlJUrtpuL2XVE7z11iTfbiGwdzJIfx24Y4amfDQmGVOktRuhg+HadPglVfyTqJOp2jCwzFv/AiiC+fww4X7a2jCQ2OWOUlSu1lqKdh3X7jkkryTqFMpmvBwcTqEOz/almu6jqLbir1rcsJDY06AkCS1q4MPzs7QjRuXLcAvtVlhwsMf2YFTOI0H2Ybe89+s2QkPjXlmTpLUrjbaCPr0gXvuyTuJOo1Zs3iBdRjJVVzFSNbhpc+2yzInSSqDgw+G3/629eOkUry7+vrszs2czsnswJ8W7qjRCQ+NWeYkSe1u5Ei4806YPTvvJKp28+fDPr1vZ9duf+RQJizcUcMTHhqzzEmS2l3v3jBsGFxxRd5JVJWKZq7+cIVLCOCXE3tnEx2c8PAFljlJUlk0DLWmlHcSVZWimasXpbHcNWcLrn6xjm5dE7z8MixYkP20yH3GMidJKotttoFPP4WHH847iapKo5mrN7M7vT98PduuJlnmJEllEeFECC2GopmrV7OvM1dLYJmTJJXNAQfAjTfC++/nnUTVwpmri84yJ0kqm5VXhh13hGuuyTuJKlbRZIf5A9ZmH6525uoissxJksrKoVY1q2iyAynxw1lHE6++yi/HPOvM1UXg7bwkSWU1eHD23+sZM2DDDfNOo4pSmOwAcBGHche78PCCLel2e+9sxqpK4pk5SVJZde0KBx0EF1+cdxJVnMKkhj+yA6dyajZzlfec7LCILHOSpLI75BC46iqYMyfvJKoo/fs3fc9VJzssEsucJKnsVl8ddtghu0RKavDuT85i97jl8zNXneywyCxzkqQO8d3vwq9/7R0halrRzNV5A9Zh7wu+zq67JA4dcIeTHdrACRCSpA6x447wySfw0EPZ3SFUYxpmrs6dSwK+N+tYunWZwS+PewsOeDnvdFXNM3OSpA4RAYcfnp2dUw0qmrl6Lj/gz2zF1Qv2otvJP8k5WPWzzEmSOszo0XDbbfCf/+SdRB2uMEN1KrvzS47hFr7Fssxx5mo7sMxJkjpM796w554wcWLeSdTh+vfncTZiDBO5kW/Tn1c+2662scxJkjrU4YfDRRfBp5/mnUQd6dVjzmFYTOX/OJzNeCzb6MzVdmGZkyR1qI03hn794NZbWzioaNYjAwe6pkk1Kvo7/G//r7L7L7fniL1ms+eAac5cbWfOZpUkdbiGZUqGDm1iZ9GsRyC7b+fYsdlz/8NfHYr+Dj+lC/Wv/IxBXW/h+J92hWtezjtdpxOphhb8qaurS9OmTcs7hiTVvI8+yi6VevhhWHvtRjsHDswKXGMDBni/zmpR9Hd4LGcznU24g13pMWA1/w4XQURMTynVtXacw6ySpA635JLwne80s0xJc7MbnfVYPQp/VxM4hKkM5XpG0IN5/h2WiWVOkpSL734XLr20ifu1Nje70VmP1aN/f+5mJ07mdG7lm6zAO59tV/uzzEmScjFgQHZXiEsvbbRj/PhslmMxZz1WlWe/ewGjuIpr2Zt1eTHb6N9h2VjmJEm5+cEP4H//t9EyJfX12SzHAQOc9VgNGs08nv1/1/Gti3bn7LEvsO2AWf4ddgAnQEiScpMSbLFFdqenJme2qrI1mnn8EUuwU5d72WH3Xvz0pvVzDlf9nAAhSap4EfD978O55+adRIul6H6rCwhGcxmrL5jF6X/dPedgtcUyJ0nK1Z57wgsvwBNP5J1Ei6xoduqPOZPXWI3LGE2XV5pYWkZlY5mTJOWqe3c48kg477y8k2iRFWan/prD+T3DuInhLMnHzlrtYLmUuYg4IyJmRMQTEXFnRKzWxDE7FPY3PD6KiOGFfZdGxD+L9g3q+G8hSWovY8fC1Knw73/nnUSLZPx4bu4xgjM4idv4BivytrNWc5DXmbmzU0obppQGAbcAJzc+IKV0b0ppUOGYHYG5wJ1FhxzXsD+l5Ml5SapiK6wAo0bB+efnnUSLYtqX6jloiSu5aZXDWStedtZqTnIpcyml94te9gRam1K7J3BbSmlu+VJJkvJ0zDFZD3jvvbyTqEmNliD557k3MXQo/PaKJdn89ZtgwYLsVl0WuQ6X2zVzETE+Il4B6mnizFwj+wJXNdo2vjBUe25ELFGWkJKkDrPmmrDrrnDxxXkn0Rc0LEEycyakxDsz32O3Y77MCbs8xrBheYdT2daZi4i7gVWa2DUupfT7ouNOAJZMKZ3SzO9ZFZgBrJZSmle07d9AD2AC8FJK6fRm3j8WGAvQv3//TWY2dfNmSVJFmDEDhgyBf/wju3+rKsTAgVmRAz6mB4O5k02YzjkDzs/OxqksSl1nLvdFgyNiAHBrSqnJ1QUj4mjgaymlsc3s3x44NqX0rdY+y0WDJanyffObMHw4HHJI3kn0mS5dICUWEOzHlXxCD65lb7oE2fCqyqKiFw2OiHWLXg4Fnmvh8JE0GmItnJkjIgIYDjzd3hklSfk4/ng466xGt/hSvgpLjZzIT3mZgVzB/nQhuQRJhcjrmrkzI+LpiJgBDAaOBoiIuoj4bcNBETEQWAO4r9H7J0XEU8BTQB/gpx0RWpJUfl//OvTtCzfc0I6/tNHF+0ya1I6/vAaMH89F3b/HFPZiKkNZio9cgqSC5D7M2pEcZpWk6jB1Kpx6Kkyfnt3yq00a3T8UyIqIS2iU7IYb4Mgxc3lg6SGs/fqD2Rm58eP98yuzih5mlSSpJbvvDinBzTe3wy8run/oZ+bOzbbrixqdxbzvxLs47DC45Z6lWfvV+12CpAJZ5iRJFScCTjklOzvX5gGkovuHlrS9ljVagmTGzGXZa/wgJo+5h403zjucmmOZkyRVpGHD2unsXHMX6Xvx/hcVncV8mQHsxh+4gCPZ+aoxOQdTSyxzkqSK1G5n58aPz66RK+bF+00rnK18kxXZlTv4EWexD9d6FrPCWeYkSRWrXc7O1ddnkx0GDMgaovcPbV7//nzA0nyTW9mDGziKCz7brsrlbFZJUkW76SY4/fR2mtmqFs27bDLDDlqRVRa8ykTGEODM3xw5m1WS1CkMG5aVuOuvzztJJ9HMmnsLFsDBfxxFlw3XZ0L/8YRnMatGt7wDSJLUkgg480w44ois2HXvnneiKtZ4zb2ZM2HsWFKC7z9az0svwR0P9qNbz5fyzalF4pk5SVLF22WX7CTRxIl5J6lyzay5d8qRb/LAA3DLLdCzZz7RtPgsc5KkqnDmmdm1c//9b95JqlgTs1J/yQ+59r1dueMO6N07h0xqM8ucJKkqbLIJbLcdnHde3kmqWKNZqb/hYH7Fkdzd70BWWimnTGozy5wkqWqMH5+Vudmz805SpYrW3LuafTiVU7lryaGs/vPv5RxMbWGZkyRVjbXWgv32g5NPzjtJlSqsuXdr3wP5Pudx+6oHsc5vf+xs1SpnmZMkVZVTToEbb4S//jXvJBWqmaVHGtyzSj3f4RKm/mUVNnjtDotcJ2CZkyRVleWXhzPOgKOOauNtvjqjhqVHZs7M/nAKS480FLp774WRI+G662CzzXLOqnZjmZMkVZ2DDspW2Lj66ryTVJhmlh5h3Djuuw/22QeuvRa23TafeCoPy5wkqep07Qrnnw/HHedSJZ/TxNIjAA/M7M9ee2Xld/vtOzaSys8yJ0mqSltvDTvsAKed1gEf1sp1aBWj0dIjAA+xFSO63MDkybDjjjlkUtlZ5iRJVeuXv4TLL4fp08v4Ia1ch1ZRipYeAXiEzfk2N3HlcU+y88455lJZWeYkSVVrpZXgF7+AMWNg3rwyfUgL16FVnMLSIwwYwP1sy9Aut3DZcU8z+Myd8k6mMrLMSZKq2n77wSqrZGfpyqKZ69Ca3V5OpQz31tdzx8UvM6LPfVx1Zx++cdYOHZ1SHcwyJ0mqahFw0UXZGbrnny/DBzRxHVqL28ulxOHem26C/ffPfu7kCbmaYJmTJFW9gQOziRD77VeG4dZG16EB2evx49v5g1pRwnDvVVfBYYfBbf+/vfsPlrK67zj+/oCxRsUEiqLyQyhQCcYJP66YFu1go0hIsKgxamyg0zEOJoy0hjGNdtQJY8bGNs6QTjQIzsRJTaNxGBxEASOQUEG5oHARaVDEij+ICrFi0lHg2z/O47jc7Lq7ZvfZu7uf18zO3T179ux3D4e93/uc55zn4bRAxNqDkzkzM2sJX/86HH98HVa3FpyHhpR+LliQ/5UTykz3LlwIc+fCo4/C+PE5xmUNp2ij7bM7Ojqis7Oz0WGYmVmd7NkDY8bAz37WghvjDh2apla7iSGncMvXdrFoEaxYASNH5h+a1YekjRHRUa6ej8yZmVnLGDAgHaGaMQP27Wt0NFWoZGFDkeneAx/vw6zhK3ngAXj8cSdy7crJnJmZtZQvfAGmT0+LAA4damAglW40XOk+dt2me383+FQuGrWNnb1HsmYNnHRSvT+Q9VSeZjUzs5bz3ntpJWduV4jo7v0ErXDBwtFHFz/XrsT0KaecArt2FW3+jTdg2jQYMQIWLYIjj6xZ5NaDeJrVzMza1sc+BvffD3ffDQ8+2IAAqtlouMp97Lq6YMKElKjec48TOXMyZ2ZmLWrAAPj5z+HKK2HLlpzfvJoErYp97JYsSddXnTcPvvvdtLjWzMmcmZm1rDPPhB/8AKZOLT6TWTfVbDRcwT52EenhN74By5blvyuK9WxO5szMrKVdeilcdx2cfz68+WZOb1rNRsNl9rHbty8t6HjwQXjySTjjjBzit6biZM7MzFreNdekhOjzn89py5JqNxq+4oq02OHQofQzq7dhA4wbB8OGwa9+BSefnEPs1nS8mtXMzNpCBFx7LaxZkzbX7d+/0RGVdugQzJ+fzou74w64+OJGR2SN4NWsZmZmBST4/vfTdOs558CrrzY6ouJ27Urbqtx3H6xb50TOynMyZ2ZmbUNKR7suvzwtjnjqqUZH9IEIuOuudE7c1KlpWnX48EZHZc3giEYHYGZmlicJrr8+Xfpq8mS4887GH/3q6oLZs+Gdd2D1ajjttMbGY83FR+bMzKwtXXIJPPIIfPObMGtWSqTy9tvfwpw5aVr1ssvgiSecyFn1nMyZmVnbGj8eNm9OF2cYOzYtjsjD/v1punfkSPj972HbNrj6aujdO5/3t9biZM7MzNraJz6RLot1660wY0Y6YvfCC/V5rz170tUbhg+HrVth7dq0Y0lPXllrPZ+TOTMzM+Cii+DZZ+H006GjA77yFdi06Y9v9+BBWLUKvvpVGDUKXnoJHnsM7r0XTj31j2/fzMmcmZlZ5iNWtogAAAnKSURBVOij4cYbYefONO164YXw6U+nCzds2AAHDlTWzp49aWuRWbNg4ECYOze19/zz6Uicz4uzWvKmwWZmZiUcOgSPP54Ss1Wr0vVdR4+GESPgxBOhT59U7+234a23UhK4fXs6B+/ss2HSJJg2LZ0bZ1atSjcN9tYkZmZmJfTqBWedlW4Ae/emqdgdO+D119NChggYMCAlbF/+cppKHTgwvdYsD07mzMzMKtSvH0ycmG5mPYX/bjAzMzNrYk7mzMzMzJpYw5M5SXMlhaSiu+xImilpR3abWVA+XlKXpOckzZek/KI2MzMz6xkamsxJGgycB/xPief7ATcBZwITgJsk9c2evgO4ChiZ3abUPWAzMzOzHqbRR+ZuB64DSu2Pcj6wMiL2RsQ+YCUwRdJJwHERsS7S3ir3ANNzidjMzMysB2lYMifpAuDliNj8IdUGAi8VPN6dlQ3M7ncvNzMzM2srdd2aRNKjwIlFnroBuB6YXK6JImXxIeXFYriKNB3LkCFDyrydmZmZWXOpazIXEecWK5d0OjAM2JytWxgEbJI0ISJeK6i6G5hU8HgQsDorH9St/JUSMSwAFkC6AsRH+RxmZmZmPVVDplkjoisiToiIoRExlJScjeuWyAEsByZL6pstfJgMLI+IV4G3JX02W8U6A1iS52cwMzMz6wkavQDiD0jqkLQQICL2AvOADdntO1kZwNXAQuA54Hng4QaEa2ZmZtZQSotB20NHR0d0dnY2OgwzMzOzsiRtjIiOcvV63JE5MzMzM6uckzkzMzOzJuZkzszMzKyJOZkzMzMza2JO5szMzMyamJM5MzMzsybmZM7MzMysibXVPnOSXgdebHAY/YE3GhxDO3F/58v9nT/3eb7c3/lq9/4+JSKOL1eprZK5nkBSZyUbAFptuL/z5f7On/s8X+7vfLm/K+NpVjMzM7Mm5mTOzMzMrIk5mcvfgkYH0Gbc3/lyf+fPfZ4v93e+3N8V8DlzZmZmZk3MR+bMzMzMmpiTuTqTdLOklyU9nd2mlqg3RdJ/S3pO0j/lHWerkTRXUkjqX+L5gwX/Jg/mHV+rqaC/Z0rakd1m5h1fq5A0T9KWbNyukHRyiXoe3zVSRZ97jNeApNskbc/6fLGkT5aot0tSV/bv0pl3nD2Np1nrTNLNwP6I+NcPqdMb+DVwHrAb2ABcHhHbcgmyxUgaDCwERgHjI+IP9iiStD8ijs09uBZUrr8l9QM6gQ4ggI1ZvX15x9rsJB0XEf+b3b8GGB0Rs4rU8/iukUr63GO8diRNBh6LiAOS/gUgIr5VpN4uoKPY93s78pG5nmEC8FxE7IyId4H/BP6mwTE1s9uB60hfqlZ/5fr7fGBlROzNfrmtBKbkFVwreT+pyByDx3jdVdjnHuM1EhErIuJA9nA9MKiR8TQLJ3P5mJ0dMr5bUt8izw8EXip4vDsrsypJugB4OSI2l6l6lKROSeslTc8jtlZUYX97fNeQpFskvQRcAdxYoprHdw1V0Oce4/Xx98DDJZ4LYIWkjZKuyjGmHumIRgfQCiQ9CpxY5KkbgDuAeaSBNw/4N9IAPayJIq/1X9wllOnv64HJFTQzJCJekfRnwGOSuiLi+VrG2Spq0N8e31X4sP6OiCURcQNwg6RvA7OBm4rU9fiuQg363GO8CuX6O6tzA3AA+I8SzUzMxvgJwEpJ2yPil/WJuOdzMlcDEXFuJfUk3QUsLfLUbmBwweNBwCs1CK0llepvSacDw4DNkiD14yZJEyLitW5tvJL93ClpNTAW8C+7ImrQ37uBSQWPBwGr6xJsC6j0+wS4F3iIIsmcx3d1atDnHuNVKNff2QKSLwKfixIn9heM8d9IWkw6XaltkzlPs9aZpJMKHl4IbC1SbQMwUtIwSUcClwFegValiOiKiBMiYmhEDCV9wY7rnshJ6ivpT7L7/YGJgBebVKnS/gaWA5Ozfu9LOpK3POdwW4KkkQUPLwC2F6nj8V1DlfQ5HuM1I2kK8C3ggoj4XYk6x0jq8/59Un8X+93aNpzM1d/3suXTW4BzgH8EkHSypGUA2cmes0n/+Z8F7ouIZxoVcCuS1CFpYfbwU0CnpM3AKuBWrxyurcL+joi9pFMMNmS372RlVr1bJW3Nvk8mA3PA47vOyva5x3hN/TvQhzR1+rSkO+Hw35nAAGBtNsafBB6KiEcaE27P4K1JzMzMzJqYj8yZmZmZNTEnc2ZmZmZNzMmcmZmZWRNzMmdmZmbWxJzMmZmZmTUxJ3Nm1tQk7a9RO6OyrRCekjS8Fm2ameXByZyZWTIdWBIRYwsvfaXE35Vm1mP5C8rMWkKWdN2WbfDaJenSrLyXpB9KekbSUknLJH2p22unAv8AXClplaShkp6V9ENgEzBY0mRJ6yRtknS/pGOz106RtF3SWknzJS3Nym+WNLfgPbZKGprd/1tJT2ZHAn8kqXdWvj+7qPtmSeslDcjKB0hanJVvlvSXkuZJmlPQ/i2SrqlfD5tZT+VkzsxaxUXAGOAzwLnAbdnl9C4ChgKnA1cCf9H9hRGxDLgTuD0izsmKTwXuiYixwDvAPwPnRsQ4oBO4VtJRwF3ANOBsil88/DCSPgVcSrpQ+BjgIHBF9vQxwPqI+AzpOpNfy8rnA2uy8nHAM8AiYGbWZi/SZQBLXZTczFrYEY0OwMysRs4CfhoRB4E9ktYAZ2Tl90fEIeA1SasqbO/FiFif3f8sMBr4L0kARwLrgFHACxGxA0DST4CryrT7OWA8sCFr6+PAb7Ln3gWWZvc3Audl9/8amAGQfb63gLckvSlpLOnyRk9FxJsVfjYzayFO5sysVajK8nLe6dbGyoi4/LCGpTFAqWsiHuDw2Y+jCtr6cUR8u8hr3osPrrF4kPLf0QuBvyMdEby7TF0za1GeZjWzVvFL4FJJvSUdD/wV6SLca4GLs3PnBgCTPkLb64GJkkYASDpa0p8D24FhBatfC5O9XaQpUSSNA4Zl5b8AviTphOy5fpJOKfP+vwCuzur3lnRcVr4YmEI6Arn8I3wuM2sBTubMrFUsBrYAm4HHgOsi4jXgAWA3sBX4EfAEaZqyYhHxOukI2E8lbSEld6Mi4v9I06oPSVoLvFjwsgeAfpKeJiViv87a2kY6/25F1tZK4KQyIcwBzpHURZp+PS1r611gFXBfNv1qZm1IHxzRNzNrTZKOjYj9kv6UdLRuYpbo1fp9JgFzI+KLtW67xPv1Iq22veT98/bMrP34nDkzawdLJX2StHBhXj0SubxJGk1aLLHYiZxZe/OROTMzM7Mm5nPmzMzMzJqYkzkzMzOzJuZkzszMzKyJOZkzMzMza2JO5szMzMyamJM5MzMzsyb2/6IHTzUcM8ZXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 加上泊松项？\n",
    "ADM_with_poisson = [i+2/mean_x_data for i in ADM]\n",
    "\n",
    "        \n",
    "        \n",
    "POW1=POW[0:int(nf)]\n",
    "\n",
    "POW_ps = [i+2/mean_x_data for i in POW1]\n",
    "PTF_ps = np.multiply(np.array(F1),np.array(POW_ps))\n",
    "F1_log = [math.log(i,10) for i in F1]\n",
    "PTF_ps_log = [math.log(i,10) for i in PTF_ps]\n",
    "\n",
    "#加上泊松项并分bin\n",
    "ADM_binned_ps = databin_20(ADM_with_poisson)\n",
    "F1_binned = databin_20(F1)\n",
    "'''ADM_times_f_b_ps = np.multiply(np.array(F1_binned),np.array(ADM_binned_ps))'''\n",
    "ADM_times_f_b_ps = np.array(F1_binned)+np.array(ADM_binned_ps)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(F1_log,PTF_ps_log,color=\"b\",linewidth=1)\n",
    "plt.scatter(F1_binned,ADM_times_f_b_ps,color=\"r\",linewidth=1)  \n",
    "plt.xlabel(\"log frequency\")\n",
    "plt.ylabel(\"log power*frequency\")\n",
    "plt.title(\"power density spectrum\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAHjCAYAAABxWSiLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcVXX9x/HXBxEEIZckcmNRE3cRKUtNMcWttDTTigxXyiV3c6FcUvqZay6l4pJmlGXlkprbL8GVdDAFLdw3tAzF3EgU+f7+OHd+DuMMc4e595577ryej8d5zD1nztz74Xpl3nzXSCkhSZKk+tYj7wIkSZLUMUObJElSARjaJEmSCsDQJkmSVACGNkmSpAIwtEmSJBWAoU2SJKkADG2SJEkFYGiTJEkqgJ55F1ANK6ywQhoyZEjeZUiSJHVo2rRpr6aUBnR0X0OGtiFDhtDU1JR3GZIkSR2KiOfLuc/uUUmSpAIwtEmSJBWAoU2SJKkADG2SJEkFYGiTJEkqAEObJElSARjaJEmSCsDQJkmSVACGNkmSpAIwtEmSJBWAoU2SJKkADG2SJEkFYGiTJEkqAEObJElSARjaJEmSCiCX0BYRZ0TEzIiYHhHXRsSy7dy3fUQ8HhFPRcSxta5TkiSpXuTV0nY7sF5KaQPgCeC41jdExBLAz4AdgHWAb0TEOjWtUpIkqU7kEtpSSrellOaXTqcCq7Rx22eAp1JKz6SU3gOuBr5cqxolSZLqST2MadsH+HMb11cGXmxxPqt0rU0RMS4imiKiafbs2RUuUSqgUaOyQ5LUEHpW64kj4g7gk218a3xK6frSPeOB+cCktp6ijWupvddLKU0EJgKMHDmy3fskSZKKqGqhLaW0zaK+HxFjgS8BW6eU2gpZs4BVW5yvArxcuQolSZKKI6/Zo9sDxwA7p5TmtnPbg8CnImJoRPQCvg7cUKsaJUmS6kleY9ouAPoDt0fEwxFxEUBErBQRNwOUJiocDNwK/AP4XUrpsZzqlSRJylXVukcXJaW0RjvXXwZ2bHF+M3BzreqSJEmqV/Uwe1SSJEkdMLRJkiQVgKFNkiSpAAxtkiRJBWBokyRJKgBDmyRJUgEY2iRJkgrA0CZJklQAhjZJkqQCMLRJkiQVgKFNkiSpAAxtkiRJBWBokyRJKgBDmyRJUgEY2iRJkgrA0CZJklQAhjZJkqQCMLRJkiQVgKFNkiSpAAxtkiRJBWBokyRJKgBDmyRJUgEY2iRJkgrA0CZJklQAhjZJkqQCMLRJkiQVgKFNkiSpAAxtkiRJBWBokyRJKgBDmyRJUgEY2iRJkgrA0CZJklQAhjZJkqQCMLRJkiQVgKFNkiSpAHrmXUARzZsH77+fPY748Hrz485eW5znkCRJ3YuhbTH88Ifw859DSh9ea37c2Wvl3t+eSofB5sc9eix8RHz02qKuL87PtLy+xBLQsycsueTCRyWu9eoFfftCnz7Z19aPl1xy0e+5JEl5MLQthtNPz468VDIMtnctJViwYOGjrWuLur44P9N8ff787Hj//Q+P1uetr737bvv3tbz+3nvw3//C3LnZ0fLx3LnZe9BWmGt+3Py1Xz9YdllYbrnsaPm4+VhmmSwsSpLUVf46KSC7S6vr/fc/GubaevzWW/Cf/8Drr8Mzz2Rfm4/m62++mQW81sFu4EBYaSVYeeXsa/Pjj3/c/66SpLYZ2qRWllwyayFbZpmuP9eCBVm4axnoXn8dXnkFXn4Z7r0XXnope/zyy/D227Diih+Gueavq64Ka66ZHf37d70uSVLxGNqkKurR48MAOGRIx/f/97/wz39mAa45zL30Ejz4IDz5ZHYss0wW3oYN+zDIDRsGQ4c6Hk+SGpmhTaojffrAaqtlR1sWLMhC3BNPwOOPZ1/vuCP7+tJLMHgwrLcefPrT8OnXRzCy/+NUoMFQklQHDG1SgfTokXWVrroqbL31wt+bNw+efhoeeSRrmTvxub15+O01WGWtUogrHcOHZ+FQklQskTpaU6KARo4cmZqamvIuQ8rXqFHMT0vw9/P/lwcf5P+Pf/wD1l4bttkGtt0WNtsMlloq72IlqfuKiGkppZEd3WdLm9TAesYHbLABbLAB7Ltvdu3dd6GpKetW/eEPYcaMLLhtu212rLuuM1glqR65jZXUzSy1FGy+OZx0Etx3H7z4Iowbl42L22mnbMbqXnvBH/+YBTxJUn0wtEnd3LLLwq67wkUXZevN3XVXNvbtgguy5Uf22gtuvTVbnFiSlB9Dm6T/FwFrrAEHHQR/+Qs89lg2ceGEE7L14g48EO6+O5vFKkmqLUObpHattBIcdhj89a9w//2wyipZcBsyBCZMgNmz865QkroPQ5uksqy+Ohx/fDZx4YYb4Nlns4V999oLpk3LuzpJanyGNkmdNnw4XHopPPVUtnzIrrtmM1Cvvjrbu1WSVHmGNkmL7eMfh2OOyRb1PeoomDgx207r/POzxX4lSZVjaJPUZT17wi67ZJMXbrgBbrsNPvWpLMTZ8iZJlWFok1RRI0bAn/4E11wDv/99tpn9FVe4ZIgkdZWhTVJVbLJJ1uJ25ZXwi19kOy384Q/QgDvnSVJNGNokVdXnPw+TJ2eL9Z54Imy/fbb7giSpcwxtkqouAkaPhr/9LQttm26aLR/yzjt5VyZJxWFok1QzSy4Jhx8O06fD88/DOuvYZSpJ5TK0Saq5lVaCSZOy8W4nnghf/jK88kreVUlSfTO0ScrNqFHw0EOw3nqw4Ybwxz/mXZEk1S9Dm6Rc9eoFP/5xFti+/30YOxbeeCPvqiSp/uQS2iLijIiYGRHTI+LaiFi2nfuei4gZEfFwRDTVuk5JtbPppvDww9C3L2ywQbZQryTpQ3m1tN0OrJdS2gB4AjhuEfdulVIanlIaWZvSJOWlXz+48EK46CLYc084+WRYsCDvqiSpPuQS2lJKt6WUmtdHnwqskkcdkurTDjtAUxPccQfstBPMmZN3RZKUv3oY07YP8Od2vpeA2yJiWkSMW9STRMS4iGiKiKbZs2dXvEhJtbXiilkX6bBhMHJktsabJHVnVQttEXFHRDzaxvHlFveMB+YDk9p5ms1SSiOAHYCDImKL9l4vpTQxpTQypTRywIABFf2zSMrHkkvC2WfDaafBtttm22FJUnfVs1pPnFLaZlHfj4ixwJeArVNqe2nNlNLLpa//johrgc8Ad1W6Vkn1bffds2VBdt0VHn0UzjgDetRDP4Ek1VBes0e3B44Bdk4pzW3nnqUjon/zY2Bb4NHaVSmpnqyzDtx/f7au29e+BnPb/JtDkhpXXv9WvQDoD9xeWs7jIoCIWCkibi7dMxC4JyIeAR4Abkop3ZJPuZLqwXLLwS23ZMuCbLWVuyhI6l6q1j26KCmlNdq5/jKwY+nxM8CGtaxLUv3r3Rt++Us46ST43Ofgpptg7bXzrkqSqi+X0CZJXRGRreG22mrZVljXXw+f/WzeVUlSdTmUV1JhjR0LV1wBO+8Md96ZdzWSVF2GNkmFtsMO8LvfwR57wM03d3y/JBWVoU1S4Y0aBX/6E+y9N/z+93lXI0nV4Zg2SQ1hk03gttuylre5c+Hb3867IkmqLEObpIax4YbZ1ldbb53tpvCNb+RdkSRVjqFNUkNZay249VYYPRp69YKvfjXviiSpMhzTJjWiSZNg6lSYMgWGDMnOu5H11oM//xkOPDAb6yZJjcDQJjWaSZNg3DiYNy87f/757LybBbfhw+HGG2HffbOxbpJUdIY2qdGMH//RjTnnzs2udzOf/jRcdx1861tw3315VyNJXWNokxrNCy907nqD23RTuOoq2GUX+Pvf865GkhafoU1qNIMGde56N7DddnDWWbD99vDii3lXI0mLx9AmNZoJE6Bv34Wv9e2bXe/GvvUtOOSQLLjNmZN3NZLUeYY2qdGMGQMTJ0Lv3tn54MHZ+Zgx+dZVB446KgttO+8M//1v3tVIUudESinvGipu5MiRqampKe8ypHyNGpV9nTw5zyrqzoIF8M1vZo9//Wvo4T9dJeUsIqallEZ2dJ9/XUnqVnr0gF/8IlsJ5aST8q5GkspnaJPU7fTpky0FctVV8Ktf5V2NJJXHbawkdUsDB2a7JWy1VbZpxOab512RJC2aLW2Suq311sta2772NXjuubyrkaRFM7RJ6ta23x6OPjpbfLf1RhKSVE8MbZK6vcMPh3XXhf33hwacUC+pQRjaJHV7EXDJJTBzJpx9dt7VSFLbnIggSWQzSq+9FjbZBDbYAEaPzrsiSVqYLW2SVDJoEPzmN7Dnnu5RKqn+GNokqYVRo7IxbrvvDu+9l3c1kvQhQ5sktXL00TBgAHz/+3lXIkkfMrRJUis9esCVV8L118M11+RdjSRlDG2S1IbllssC24EHwpNP5l2NJBnaJKldI0fCiSfCN77h+DZJ+TO0SdIiHHQQrLgi/PCHeVciqbsztEnSIkTA5ZfDr34Fd9yRdzWSujNDmyR1YMAAuOIK2GsvePXVvKuR1F0Z2iSpDKNHZ2Pb9t3X/Ukl5cPQJkllmjABZs2Ciy7q4hONGpUdktQJ7j0qSWXq1Qt+/WvYfHPYYgtYd928K5LUndjSJkmdMGwYnHZa1lX67rt5VyOpOzG0SVIn7bNPFt7c5kpSLRnaJKmTImDiRLjuOrj99ryrkdRdGNokaTEstxxccgnstx+88Ube1UjqDgxtkrSYttsOtt8ejjwy70okdQeGNknqgjPPzHZK+POf865EUqMztElSF/Tvn21zNW4cvP563tVIamSGNknqoi98Ab78ZTjssLwrkdTIDG2SVAGnnQb33gs33JB3JZIalaFNkiqgXz/4xS/ggAPgtdfyrkZSIzK0SVKFfP7zsPvu8L3v5V2JpEZkaJOkCpowAZqasoV3JamSDG2SVEF9+2aL7h58sIvuSqosQ5skVdiWW8IXvwjHHpt3JZIaiaFNkqrgJz+BP/0J7rkn70okNQpDmyRVwbLLwnnnwf77w7vv5l2NpEZgaJOkKtl1V1h7bfjxj/OuRFIjMLRJUhVdcAFcdBE8+mjelUgqOkObJFXRSivBqafCfvvBBx/kXY2kIjO0SVKV7bcf9OoFP/953pVIKrKeeRcgSY2uR49s7bbNNoNddoFV8i5IUiHZ0iZJNTBsGBx0EBxxRN6VSCoqQ5sk1cixx8K0aXDrnE/nXYrUeEaNyo4GZmiTpBrp0yebTXrQk4fx7oJeeZcjdV8FDXiGNkmqoR12gOH9nuInL3wj71IkFYyhTZJq7JzVL+D8l3blqafyrkRSkeQW2iLilIiYHhEPR8RtEbFSO/eNjYgnS8fYWtcpSZW26lKzOXbQrzn4YEgp72okFUWeLW1npJQ2SCkNB24ETmh9Q0QsD5wIbAJ8BjgxIparbZmSVHmHrvx7XnoJ/vCHvCuRVBS5hbaU0pstTpcG2vr35nbA7SmlOSml14Hbge1rUZ8kVdOSPT7gwgvh8MPhrbfyrkZSEeQ6pi0iJkTEi8AY2mhpA1YGXmxxPqt0ra3nGhcRTRHRNHv27MoXK0kVtvnmsM02cNJJeVciFUBBZ3xWUlVDW0TcERGPtnF8GSClND6ltCowCTi4rado41qbI0BSShNTSiNTSiMHDBhQuT+EJFXR6afDVVfB9Ol5VyKp3oNhVbexSiltU+atvwZuIhu/1tIsYFSL81WAyV0uTJLqxIABcMopcMABcPfd2ZZXktSWPGePfqrF6c7AzDZuuxXYNiKWK01A2LZ0TZIaxv77w/z5cMUVeVci5ajOW7nqQZ7/pjut1FU6nSyMHQoQESMj4lKAlNIc4BTgwdLxo9I1SWoYPXrAhRfCccfB66/nXY2kelXV7tFFSSl9tZ3rTcB+Lc4vBy6vVV2SlIcRI2DXXeHEE+G88/KuRurmmlv8Jk/Os4qP6LClrbRWmiSpyk49Fa6+GmbMyLsS7KpSdfn5WizldI/+NSKuiYgdI6Kt2ZySpAr4+Mez5T8OOcSdEiR9VDmhbU1gIrAn8FRE/Dgi1qxuWZLUPY0bB3PmwDXX5F2JpHrTYWhLmdtTSt8gG2s2FnggIqZExOeqXqEkdSM9e2Zj2o46Ct55J+9qJNWTcsa0fTwiDo2IJuAo4HvACsCRZOurSZIqaMstYbPN4Cc/ybsSKUeOe/uIcrpH7wc+BnwlpfTFlNIfU0rzS7M8L6pueZLUPZ1+OvzsZ/DMM3lXIqlelBPahqWUTkkpzWr9jZSS/w6UpCpYdVU44gg48si8K5FUL8oJbbdFxLLNJ6XdCdyVQJKq7Mgjsz1Jb7st70ok1YNyQtuAlNJ/mk9SSq8Dn6heSZIkgKWWgnPOyZYAee+9vKuRlLdyQtsHETGo+SQiBgOuICRJNbDTTjBkCFxwQd6VSMpbOdtYjQfuiYgppfMtgHHVK0mS1CwCzj03m006ZgwMHJh3RZLyUs46bbcAI4DfAr8DNk4pOaZNkhbHpEkwdSpMmZI1oU2a1OGPDBsG3/42/PCH1S9PWojLbtSVcrpHAXoDc4A3gHUiYovqlSRJDWrSpGzLg3nzsvPnn8/OywhuJ5wA118PDz9c5RqlRtEycDZI+OywezQifgLsATwGLChdTsBdVaxLkhrP+PEwd+7C1+bOza6PGbPIH1122Wxf0sMOgzvvzLpNJXUv5bS0fYVsrbYvppR2Kh07V7swSWo4L7zQueut7L8/vPYaXHddBWuSVBjlhLZngCWrXYgkNbxBgzp3vZWePeHss7N9SZt7WCV1H+WEtrnAwxFxcUSc13xUuzBJajgTJkDfvgtf69s3u16m0aNh3XWzTeUldS/lLPlxQ+mQJHVF87i1fffNmsoGD84CWwfj2Vo780zYdNNsRqlLgEjdR4ehLaV0ZUT0AQallB6vQU2S1LjGjIFLLskeT568WE+x5pofLgEycWLlSpNU3zrsHo2InYCHgVtK58MjwpY3ScqRS4BI3U85Y9pOAj4D/AcgpfQwMLSKNUmSOtC8BMjhh0NyY0GpWygntM1PKb3R6pp/RUhSzvbfH2bPzlrcJDW+ckLboxHxTWCJiPhURJwP3FfluiRJHejZE845xyVApO6inND2PWBdYB7wG+BN4LBqFiVJKs/o0bDOOi4BInUH5cwenQuMLx2SpDrjEiCqa817fi7mbGl9qJy9R++kjTFsKaUvVKUiSZXhX5DdRvMSICecABdfnHc1kqqlnMV1j2rxeCngq8D86pQjSVocJ5wAw4bBgQfChhvmXY2kaiine3Raq0v3RsSUKtUjSVoMzUuAHHYY/OUvEJF3RZIqrZzFdZdvcawQEdsBn6xBbZKkTnAJEKmxldM9Oo1sTFuQdYs+C+xbzaIkSZ3XvATIAQfADjtA7955VySpksrpHnX3A0kqiNGjYe214fzzs/XbpJqq5UzRbjgrtZzZo7su6vsppT9WrhxJUledeSZsvjmMHQsDBuRdjaRKKad7dF9gU+AvpfOtgMnAG2TdpoY2Saojw4bBmDHZjNILL8y7GqnOFajFrpzQloB1Ukr/BIiIFYGfpZT2rmplkqTFdsIJsNZa2RIg66+fdzWSKqGcbayGNAe2kleANatUjySpApZfPgtuRxwB6SPLo0tdNGrUhy1UqplyQtvkiLg1IvaKiLHATcCdVa5LktRF3/kOzJoFN92UdyWSKqHD0JZSOhi4CNgQGA5MTCl9r9qFSZK6Zskl4eyz4cgj4b338q5GUleV09IG8BBwU0rpcODWiOhfxZokSRWyww4wdKgTEqRGUM6OCPsDvweatyFeGbiumkVJkirnrLPg1FPhtdfyrkRSV5TT0nYQsBnwJkBK6UngE9UsSpJUOeuuC3vske1NKqm4yglt81JK/z8aIiJ6ki0DIkkqiJNOgquvhn/8I+9K2uFsRKlD5YS2KRFxPNAnIkYD1wB/qm5ZkqRKWmEFOP74bFKCpGIqJ7QdC8wGZgDfAW4GflDNoiRJlXfQQfDUU3DLLXlXok6xFVIli9wRISKWAK5MKX0LuKQ2JUmSqqFXr2xf0iOOgG22gZ7l7IkjqW4ssqUtpfQBMCAietWoHklSFe20E6y0Elx8ccf3Sqov5fw76zng3oi4AXin+WJK6exqFSVJqo6IbMHdbbaBb34Tllsu74oklaucMW0vAzeW7u3f4pAkFdAGG8Auu8App+RdiaTOaLelLSKuSintCfwnpXRuDWuSJFXZKafAOuvAd78La66ZdzVSByZNgqlTYd48GDIEJkzIu6JcLKqlbeOIGAzsExHLRcTyLY9aFShJqrxPfAK+/3046qi8K5E6MGkSjBuXBTaA55/Pzl95pTqv98orWUCcMiULiJMmVed1FsOiQttFwC3AWmR7j05rcTRVvzRJUjUdeig89hjccUfelXRTLuVRnvHjYe7cha/NnQvPPlv513rlFXjiiY8GxDoJbu2GtpTSeSmltYHLU0pDWx2r1bBGSVIV9O4Np58Ohx8O8+fnXY26pJED4AsvtH29OVhV0rPPwoIFC1+bOzcLjnVgkRMRIuKbKaUDIuLrtSpIklQ7u+4Kyy8Pl12WdyVSOwYNavt6796Vf632gmB7wbHGOpo9unJE7A6sUotiJEm1FQHnnAMnnghvvJF3NVIbJkyAvn0Xvta3Lwwd+uF580SFro5Day8Ithcca6zd0BYRJwLLA78Glo+IE2pWlSSpZkaMgB137LYT8lTvxoyBiRM/DFSDB2fnAwdm55WcqDB0KPTIotE79OVN+mcBsU7+51jUmLaTgTnAt4A5KaUf1awqSVJNTZiQdZE+/XTelUhtGDMGPvtZ2HJLeO657LxZJScqDBwIa67JvUuOYm3+wbUf3z8LiC1fL0cddY++nFK6GnipFsVIkvKx4opw5JHZMiBSoVRwokJKcOZ7h7Arf+Dn613I2FfPqpvABh1vY/WbiHg0pbReTaqRJOXm8MNh7bVh8mQYlXcxUrkGDcq6RFvr5ESFOXNgr8cm8O/3luOBEQcweKkqrQPXBR1tGL8AeCQi6mMEniSpavr0gZ/8JAtvH6RydjlUQyj6ciHlTFTowNQ312HECPhUn1ncNfyQugxsUN7eoysCj0XE/0bEDc1HtQuTJNXe7rvD0kvDZf/cMe9SpPJ0NFFhEVKCc2btxpcfPZVzz4WzVr+QXj3qd9HCjrpHAU6uehWSpLoQAeefD9tvsg+7DZiCexaqJrq6t+iYMXDJJdnjyZOzr83n7Xj9/X7svQu8/O+tmbrRgQz98tVwTqcrr6kOW9pSSlOA54AlS48fJNvWSpLUgDbaCHZd4W5OeG6fvEtRd1DrvUWBB98cxoiHLmHwYLhn+PcY2udfC6/1NnVqVV9/cXUY2iJif+D3wMWlSysD13XlRSPilIiYHhEPR8RtEbFSO/d9ULrnYbtkJal2Th16Gb+bPYpHHsm7EjW8Gu4tmhKcdx588dHTOHO1Czn3XLLu0FdeWTg4zpuX7UFaJ3uONitnTNtBwGbAmwAppSeBT3Txdc9IKW2QUhoO3Ai0t3Dvf1NKw0vHzl18TUlSmT6+5JucPOQXfO972S86qWoWtWRHV3c4aOE/8/ux225w5ZVw/0YH8dUBd334zWef/WhwXLCgbvYcbVZOaJuXUnqv+SQiegJd+l84pfRmi9Olu/p8kqTKG7fijbz1Flx9dd6VqKF1tEVUBbpLp721JhtPu5gVV4T77oPV+7y88A11vudos3JC25SIOB7oExGjgWuAP3X1hSNiQkS8CIyh/Za2pSKiKSKmRsRXOni+caV7m2bPnt3V8iSp21siFnD++XD00fD223lXoy6p52U92lqyo7XF7C5NCS64ALafcTr/M/QSLrigneXb6nzP0WblhLZjgdnADOA7wM3ADzr6oYi4IyIebeP4MkBKaXxKaVVgEnBwO08zKKU0Evgm8NOIWL2910spTUwpjUwpjRwwYEAZfyxJysnkyR/OcKtzm2+e/a6vk60X1YhaL9nRnk7ucPCf+f342t9P5vLLs+7Q3T8xuf2bhw79aHDs0aPuPvjlhLZRwKSU0tdSSrullC5JqeMRDimlbVJK67VxXN/q1l8DX23nOV4ufX0GmAxsVEa9kqQKOv30bPWEJ5/MuxI1rJZ7iw4e3PY9ndjhoKkJRkybyCd7zeG++2CNPh3sxjlw4MLBsXdvWHPNutrCCsoLbXsBD0fE/RFxekTsFBHLdeVFI+JTLU53Bma2cc9yEdG79HgFsskQf+/K60qSOm+llbI9SQ87LO9K1C10YYeDlODcc2GHHeD01S7mgk+dy1JLlfm6LYPjZz9b1uK8tVbOOm3fTimtSdYaNgv4GVl3aVecVuoqnQ5sCxwKEBEjI+LS0j1rA00R8QhwJ3BaSsnQJkk5OOwweOopuPHGvCtRw1vMHQ5ef78fuz52ClddlS2zttuAKTUotrY63BEhIr4FfB5YH3gVuAC4uysvmlJqrzu0Cdiv9Pi+0mtKknLWq1fWgnHwwbDNNpTfeiEtjk7ucPDAA7DHQ5ew08fv4+p7O71XfGGU0z36U2A4cAlwSErp9JTS/dUtS5JUb7bfHtZdF84+O+9KylDPsyVVMSnBOefAl74EZ632c85b4/yGDWxQRktbSmmFiFgX2AKYUBqP9nhKac+qVydJqivnnAOf/jR885vZmqdSXua835+9Hz+Gf/4X/vpXGLp3lzoBC6Gcbaw+BgwCBgNDgGWABdUtS5JUj1ZbDQ4/HA45JO9KqsyWuro2dSqMmHYJqy31T+65p6w5Cg2hnO7Re4CdgOnAHimlYSmlsdUtS5JUr44+OtuW8frWCzhJVZYSnHUW7Lwz/HSNCzhnjZ/Rq1feVdVOOd2jGwBERH/cbkqSur3eveHCC2HsWNh6a+jXL++K1B289v7H2Gvmscyel008GLLXPXmXVHPldI+uFxF/Ax4F/h4R0yJiveqXJkmqV1ttlS1ndfLJeVdSR+xSrZr77ssWyx3W90Xuuqv7jqcsp3t0InBESmlwSmkQcGTpmiSpGzvzTLjiCpg+Pe9K1KgWpOCMM2CXXeD8Nc7jzNUv7Fbdoa2VE9qWTind2XySUpoMLF21iiRJhTBwIJx6Knz3u7DA6WmqsFffX4adHv0xf/xj1h268woXvHgiAAAb+0lEQVT35V1S7soJbc9ExA8jYkjp+AHwbLULkyTVv/33zwaHX3ZZ3pWoMCZNyqZ/TpmS9XNOmvSRW+65J+sOXbfvc9x1V/vbkXY35YS2fYABwB+Ba0uP965mUZKkHJXxS7VZjx5w0UUwfjzM7uoGh2p8kybBuHEwb152/vzz2XnpM7YgBaedBrvtBhd+6hxOX/1illxyMV+n5Wf4lVcq9kfIUzmzR18HDomIZYAFKaW3ql+WJCkX7f1ShWxroTZsuCHsuScccQRcdVWN6lQxjR8Pc+cufG3uXBg/ntkrbci3Zx7Pm/PhwQdh1T2nLt5rtPUZ7lFOG1X9K2f26KcjYgbwCDAjIh6JiI2rX5okqeYW8Ut1UU4+OevS+vOfq1ibiu+FF9q8fOfzq7HRtEvZsN/TTJ4Mq67ahddo6zO8YAE8W/yRXeVEz8uAA1NKQ1JKQ4CDgF9UtSpJUj7a+aXa7vWSfv1g4sRsUsJb9seoPYMGLXSagEvYjzFL/IbLh/2E01abuHjdoS2191ltbnlrqWDdqOWEtrdSSv+/oVdK6R7A/yUlqRG1+qXa4fUWRo/OFts97rgK16TGMWEC9O0LwH9YhgP4Odf0+DoPnXsP2y7fVJnXaO+z2non+Vde+Wg36hNP1HVwKye0PRARF0fEqIjYMiJ+DkyOiBERMaLaBUqSaqjFL9X/17dvdr0MZ50F114Ldzf+3t1aHGPGwMSJXNdzN4bxOEOXfYNbrvgXnzzoq5V7jbY+wz16fLhBaXPr2syZhetG7XAiAjC89PXEVtc3JWvZ/EJFK5Ik5ad5ssG++2YtEIMHZ78E25mE0Npyy8EFF2Q//sgj0KdPFWtV4bz7Lhw9dQw39vwC1613Mp/7288r/yJtfYb79MkWFmw9SaEti/pezsqZPbpVLQqRJNWJMWPgkkuyx5Mnd/rHd9kFfvMb+NGP4H/+p7Klqbgefxz22APWWAP+tvH+LNvz7eq9WOvPcPP2Ym1NUmitdTdqHWmMObCSpLpy/vlw+eXw0EN5V6K8pQRX/ms7Nt8cDjgArrmG6ga2RelgQs1C3ah1qJzuUUmSOmXgQDjjDNhnn2wLou68X2R39tZbcODM45n29pr8ZSqsv37OBQ0alE04aEvLbtQ6ZUubJKkq9twzW2/rRz/KuxLlYdo0GDEClurxHk0jvpN/YIP2JymstRY891xdBzYoo6UtInZt4/IbwIyU0r8rX5IkqRFEwKWXZjsm7LgjbLpp3hWpFlKCc3+a5aPzz4evX3Rm3iV9aFGTFAqgnO7RfYHPAXeWzkcBU4E1I+JHKSU3LZEktWngQLjwQvj2t+Hhh7NFeNW4Xn1/GfaeeQyv/Bf++ldYbTXgoryraqW9SQoFUE736AJg7ZTSV1NKXwXWAeYBmwDHVLM4SVLx7bILfP7zcOSReVeiarr9dtiw6VLW6vsC99xTCmyqqHJC25CUUsvlgf8NrJlSmgO8X52yJEmN5Nxz4bbb4Kab8q5ElTZvHhx1FOy9N1y51mmcsfpF9TXxpPVWVZMm5V3RYisntN0dETdGxNiIGAvcANwVEUsD/6lueZKkRvCxj8GVV8L++8Ps2XlXo0qZORM++1l46qms+3ub5ablXdLC2tqqaty4ut6qalHKCW3NG8QPBzYCrgQOSim948K7kqRybbFFNpzoO9/JBquruFKCi1/eic03h+9+N9u6bIUV8q6qDc8++9HFdOfOreutqhalnB0RUkTcA7xHtm3VAyn5v5skqfNOPRU22QQmTszCm4rn1Vdhv8dO5fl5A7n7AVh77dI3mrsh583LuiE7sf1Z1bS3JVUdb1W1KB22tEXE7sADwG7A7sBfI2K3ahcmSWo8vXvDb38LP/gBTJ+edzXqrDvugOHDYY0+LzF1owMXDmxtdUPmPX6svS2p6nirqkUpp3t0PPDplNLYlNK3gc8AP6xuWZKkRjVsGJx9drYP5Tvv5F2NyjFvHhx9NOy1F/ziF3Dm6hfSu0eLuYht7ek5d252PU9Dh350Md2+fet6q6pFKSe09Wi1iO5rZf6cJElt2nPPbAD7wQfnXYk6MnPuID73OXjiiWyywejRbdzU3p6eHe31WW0DB2Z98c0ta4MHZ+cFWUy3tXLC1y0RcWtE7BURewE3ATdXtyxJUqO74IJsCNQvf5l3JWrLggVw7qyvsvnfzmfcOLjuutJkg7aW0Bg0qO0nae96LY0Zk/0LYcsts62q8h5n1wUdhraU0tHARGADYENgYkrJRXUlSV2y9NJwzTXZorsPv71G3uWohRdfhG23hav//QXu3+hAvvvdbFuydseu7bhj292QEybUvPZGVlY3Z0rpDymlI1JKh6eUrq12UZKk7mG99eBnX7+bXaeN57UpMwq/+GnRpQRXXQUbbwxf+ALcvdEhfKrvSx/e0N7YtZtvbrsbssCtWvWo3SU/IuItsiU+PvItspVAPla1qiRJ3cOkSex++Tj+zd5cz86Mff5Klhg3Lvuev/Br6tVXszXXZs6EW2+FjTYCbvtg4ZsWNXat9Z6eRVan9bfb0pZS6p9S+lgbR38DmySpIkotNwfxM8bySz6gZ33MOuxmbnztc2ywQTapsqmpFNjaUi9j1xpoa6rOcBaoJCk/pZabAJZgAT1YwMNsmP+sw27irbdg/8eP4ntPHcLVV8MZZ8BSSy3iByZMyH/sWnvj6gq6NVVnGNokSflp1ULTgwUczwTu/cQuORXUfdx+O6y/PiwgeGTjfdliizJ+aMyY/MauTZ6cHe2Nqyvo1lSdYWiTJOWnVctNDxKH9p7Ibu/+imeeybGuBvbGG7D//rDffnDRRXDZsDP4WM+5Hf9gs7yX0GivFbagW1N1hqFNkpSfNlputrtsd074nz588Yvw+uv5lpe7Co/duvnmbMbuEkvAjBmw/fa1e+2KaW/8XEG3puqMDjeMlySpqtqYdXgA8PTTsNNOcNttHx1G1S20N3YLOt26Nef9/hz+9EHcfTBccQVsvXXtXrviJkzIamnZRdq3b30s5FtltrRJkurS6afD6qvD174G77/f8f0Np0L7eV53HazfdDnL9HyH6dPLCGwVfO2qaG9cXUG3puoMQ5skqS716AGXXpp15e21V7atUrfSxf08X3wRdtkFjjkGfrP2KZy3xvn061eb1666vMfV5cTQJkmqW0suCb/9LcyaBQcc0M2C22KuiTZ/Pvz0p9laa8OHw/TpsMWy02vy2h9Rr+PiCsrQJkmqa336wI03wmOPwXe+042C22KsidbUBJ/5DNxwA9x7L5x44mKOz6/EemztjYszuC02Q5skqe717w9//jM8/ni2XEXhg1s5LVDlrok2aRJv3P93Dp2yC1/aZDaHb3If//u/MGxYF+qrxHps9TwurqAMbZKkQujfP1uy4umnYc894b338q5oMXWmBaqDsVsLrprEZfvcw5j3LmNLpvD4gtXZ85ejiV9XoDWrq+PG6n1cXAEZ2iRJhdGvX9bi9s478KUvZdswFU6FWqDuuw8+s9/6/PW94dzAl9mVa1mGt+qnNate9iltIIY2SVKh9OkDv/99trn5qFHwr3/lXVEndbEF6qWX4Fvfgt13h8PfO52L+S49SIv1XFVVD/uUNhhDmySpcHr2zLZg+spXsoH3Dz6Yd0WdsJgtUP/5T9aAtsEG2RCzmTNhzOB7iMV4rprIc5/SBmVokyQVUgT88Idw3nmw445w5ZV5V1SmTrZAvbugF2e+uAdrrgn//Cf87W/Zrf36df65aq6brqdWLYY2SVKhfeUr2STMCROyTdDffjvvijpQZgvU++/D5ZfDmg9cxb1vrMfkydn5Qo1oXW3Nch21QnHvUUlS4a2zDkybBt/7HowYAb/uP4yR/R/Pu6z2tbHfarN587L9QU87LRu399u1T+Zzy/wd1pnc+lk6fK5Fquf9RdUmW9okSQ2hf/8s7JxyCuw44zSOeWbcRyZp1rN33oFzz832W73hhixT/eUvZIGtGhptHbXJkzsXWgvI0CZJaih77AHTR+7LC+8OZL314E9/gpQ6/rm8PPccHH101rM5ZQpcfz3cdBNsummVX9h11ArH0CZJajif7DWH36xzChdeCMcemy0NMnVq3lV96IMP4PY5G7PLo6ew8cZZqHzgAfjjH2HjjWtUhOuoFY6hTZLUsLbbDh55BMaOha99DbbdFu64o4OWtyoOzp85E447LmtVO/bZcWy7/IM8/zyceSastlrFXqY89T7zVB9haJMkNbSePWGffbLtr775TTj00Gyts7PPhn//u9XNFd7kPKVsiY4TT4Thw2GrrWD+fLjlFpi28Xc4YKUbsqU78lDv66h1gzFqneXsUUlSt9CrF+y1F3z723D33dnyGSefDBtuCF/8YrbW27rH/4Ae7Q3OLyPMpARPPJ410jUfSy2VLUty/vnZOLUllqjOn2+xLO7MU+XC0CZJ6lZ69MjWet1yS/jvf+HOO7OB/1/5Crz6wsOcyVHsw2UsQWIB8E9W5J3ne9N3Vha45s3Lfm727GxLqZdegqeeghl/O59H3xnKMqOz595qKzjpJFhjjWwhYHXA0NghQ5skqdvq0ydrYdtxx+z81VVH8dSsXixgCWAB79OT4/gf7u+5BXM3ybo2+/TJWs9WWAFWXhlWWgnWXRe+PvQy1lv6WVa49/pc/0xqXIY2SZJKVjjtKFYYNw7mzgdgCd7jl30PLI31GrroH/7DwzWoUN1Z7hMRIuKoiEgRsUI73x8bEU+WjrG1rk+S1I3U++B8VVedT37ItaUtIlYFRgNtruQXEcsDJwIjgQRMi4gbUkqv165KSVK34uB81am8W9rOAb5PFsjash1we0ppTimo3Q5sX6viJEmS6kVuoS0idgZeSik9sojbVgZebHE+q3StrecbFxFNEdE0e/bsClYqSZKUv6p2j0bEHcAn2/jWeOB4YNuOnqKNa222yqWUJgITAUaOHFnHu8xJkiR1XlVb2lJK26SU1mt9AM8AQ4FHIuI5YBXgoYhoHfBmAau2OF8FeLmaNUuSCq6K21BJecplIkJKaQbwiebzUnAbmVJ6tdWttwI/jojlSufbAsfVpEhJUvG0tw0VOANUbSvQZJO8JyJ8RESMjIhLAVJKc4BTgAdLx49K1yRJ+qjx47Ntp1pq3oZKKri6WFw3pTSkxeMmYL8W55cDl+dQliSpaF5ocwWp9q9LBVJ3LW2SJC22QYM6d10qEEObJKlxTJgAffsufK1v3+y6VHCGNklS43AbKjWwuhjTJklSxbgNlRqULW2SJEkFYGiTJEkqAEObJElSARjaJEmSCsDQJkmSVACGNkmSpAIwtEmSJBWAoU2SJKkAXFxXkiQ1ngZcWNmWNkmSpAIwtEmS1FWTJsHUqTBlCgwZkp3Xi3quTZ1iaJMkqSsmTYJx42DevOz8+eez83oIR/VcmzrN0CZJUleMHw9z5y58be7c7Hre6rk2dZqhTZKkrnjhhc5dr6V6rk2dZmiTJKkrBg3q3PVaqufa1GmGNkmSumLCBOjbd+Frfftm1/NWz7Wp0wxtkiR1xZgxMHEi9O6dnQ8enJ2PGZNvXVDftanTXFxXkqSuGjMGLrkke1xvi7rWc23qFEObJEmqLcPjYrF7VJIkqQBsaZMkfZQtIVLdsaVNkiSpAAxtkiRJBWBokyRJKgBDmyRJUgE4EUGSlD8nPjQG/ztWlS1tkiRJBWBokyRJKgC7RyVJUufYDZoLW9okSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQXgkh+SJNUzl9dQiaFNkiQ1tgYJvnaPSpIkFYChTZIkqQDsHpUkScXTIF2enWFokySptaIFgqLVq8ViaJMkSd1LQUOuoU2SpDwUNDgoP05EkCRJKgBDmyRJUgEY2iRJkgrAMW2SJHVnjq0rDFvaJEmSCsDQJkmSVACGNkmSpAIwtEmSJBWAoU2SJKkADG2SJEkFYGiTJEkqAEObJElSARjaJEmSCsDQJkmSVAC5hraIOCoiUkSs0M73P4iIh0vHDbWuT5IkqV7ktvdoRKwKjAZeWMRt/00pDa9RSZIkSXUrz5a2c4DvAynHGiRJkgohl9AWETsDL6WUHung1qUioikipkbEV2pRmyRJUj2qWvdoRNwBfLKNb40Hjge2LeNpBqWUXo6I1YC/RMSMlNLT7bzeOGAcwKBBgxazakmSpPpUtdCWUtqmresRsT4wFHgkIgBWAR6KiM+klP7V6jleLn19JiImAxsBbYa2lNJEYCLAyJEj7XKVJKkeTJ6cdwUNo+bdoymlGSmlT6SUhqSUhgCzgBGtA1tELBcRvUuPVwA2A/5e63olSZLqQW6zR9sSESOB76aU9gPWBi6OiAVk4fK0lJKhTZLUMVt31IByD22l1rbmx03AfqXH9wHr51SWJElSXck9tEmSpDpli2VdcRsrSZKkAjC0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQB65l2AJEkNYfLkvCtQg7OlTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQXgNlaSJDU6t9hqCLa0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNokSZIKwNAmSZJUAIY2SZKkAjC0SZIkFYChTZIkqQAMbZIkSQUQKaW8a6i4iJgNPJ9zGSsAr+ZcQ3fi+117vue15ftde77ntdWd3+/BKaUBHd3UkKGtHkREU0ppZN51dBe+37Xne15bvt+153teW77fHbN7VJIkqQAMbZIkSQVgaKueiXkX0M34ftee73lt+X7Xnu95bfl+d8AxbZIkSQVgS5skSVIBGNokSZIKwNBWIRFxUkS8FBEPl44d27lv+4h4PCKeiohja11no4mIoyIiRcQK7Xz/gxb/TW6odX2NqIz3fGxEPFk6xta6vkYREadExPTSZ/e2iFipnfv8jFdIJ95zP+MVEBFnRMTM0nt+bUQs2859z0XEjNJ/l6Za11lPHNNWIRFxEvB2SunMRdyzBPAEMBqYBTwIfCOl9PeaFNlgImJV4FJgLWDjlNJHFmWMiLdTSv1qXlyD6ug9j4jlgSZgJJCAaaX7Xq91rUUXER9LKb1ZenwIsE5K6btt3OdnvELKec/9jFdORGwL/CWlND8ifgKQUjqmjfueA0a29Xd8d2NLW219BngqpfRMSuk94GrgyznXVGTnAN8n+4tTtdHRe74dcHtKaU7pl9jtwPa1Kq6RNIeHkqXxc151Zb7nfsYrJKV0W0ppful0KrBKnvUUgaGtsg4uNfNeHhHLtfH9lYEXW5zPKl1TJ0XEzsBLKaVHOrh1qYhoioipEfGVWtTWqMp8z/2MV1BETIiIF4ExwAnt3OZnvILKeM/9jFfHPsCf2/leAm6LiGkRMa6GNdWdnnkXUCQRcQfwyTa+NR64EDiF7MN1CnAW2Ydwoado42f913M7Oni/jwe2LeNpBqWUXo6I1YC/RMSMlNLTlayzkVTgPfcz3gmLer9TStenlMYD4yPiOOBg4MQ27vUz3gkVeM/9jHdCR+936Z7xwHxgUjtPs1npM/4J4PaImJlSuqs6Fdc3Q1snpJS2Kee+iLgEuLGNb80CVm1xvgrwcgVKa0jtvd8RsT4wFHgkIiB7Hx+KiM+klP7V6jleLn19JiImAxsB/kJrRwXe81nAqBbnqwCTq1JsAyj37xTg18BNtBHa/Ix3TgXecz/jndDR+12ayPElYOvUziD7Fp/xf0fEtWRDjbplaLN7tEIiYsUWp7sAj7Zx24PApyJiaET0Ar4OONurk1JKM1JKn0gpDUkpDSH7S3RE68AWEctFRO/S4xWAzQAnfSyGct9z4FZg29J7vxxZy9ytNS63IUTEp1qc7gzMbOMeP+MVVM57jp/xiomI7YFjgJ1TSnPbuWfpiOjf/Jjs/W7r92u3YGirnNNLU5KnA1sBhwNExEoRcTNAacDlwWT/g/8D+F1K6bG8Cm5EETEyIi4tna4NNEXEI8CdwGnO1K28lu95SmkO2fCAB0vHj0rX1HmnRcSjpb9TtgUOBT/jVdbhe+5nvKIuAPqTdXk+HBEXwcK/N4GBwD2lz/gDwE0ppVvyKTd/LvkhSZJUALa0SZIkFYChTZIkqQAMbZIkSQVgaJMkSSoAQ5skSVIBGNok1b2IeLtCz7NWaWmBv0XE6pV4TkmqFUObpO7kK8D1KaWNWm71FBn/PpRU1/xLSlJhlMLVGaUFUGdExB6l6z0i4ucR8VhE3BgRN0fEbq1+dkfgMGC/iLgzIoZExD8i4ufAQ8CqEbFtRNwfEQ9FxDUR0a/0s9tHxMyIuCcizouIG0vXT4qIo1q8xqMRMaT0+FsR8UCpZe/iiFiidP3t0qbkj5Q2eR9Yuj4wIq4tXX8kIjaNiFMi4tAWzz8hIg6p3jssqZ4Z2iQVya7AcGBDYBvgjNIWcrsCQ4D1gf2Az7X+wZTSzcBFwDkppa1Kl4cBv0wpbQS8A/wA2CalNAJoAo6IiKWAS4CdgM/T9ubXC4mItYE9yDa6Hg58AIwpfXtpYGpKaUOy/RP3L10/D5hSuj4CeAy4DBhbes4eZFvftbeptqQG54bxkopkc+A3KaUPgFciYgrw6dL1a1JKC4B/RcSdZT7f8ymlqaXHnwXWAe6NCIBewP3AWsCzKaUnASLiV8C4Dp53a2Bj4MHSc/UB/l363nvAjaXH04DRpcdfAL4NUPrzvQG8ERGvRcRGZNv5/C2l9FqZfzZJDcbQJqlIopPXO/JOq+e4PaX0jYWeOGI40N5+f/NZuMdiqRbPdWVK6bg2fub99OH+gR/Q8d/DlwJ7kbXwXd7BvZIamN2jkorkLmCPiFgiIgYAW5BtIn0P8NXS2LaBwKjFeO6pwGYRsQZARPSNiDWBmcDQFrNNW4a658i6MomIEcDQ0vX/BXaLiE+Uvrd8RAzu4PX/FzigdP8SEfGx0vVrge3JWhRvXYw/l6QGYWiTVCTXAtOBR4C/AN9PKf0L+AMwC3gUuBj4K1n3YtlSSrPJWrR+ExHTyULcWimld8m6Q2+KiHuA51v82B+A5SPiYbLA9UTpuf5ONj7uttJz3Q6s2EEJhwJbRcQMsm7TdUvP9R5wJ/C7UreppG4qPmyll6Tiioh+KaW3I+LjZK1vm5UCXaVfZxRwVErpS5V+7nZerwfZ7NavNY+rk9Q9OaZNUqO4MSKWJZtAcEo1AlutRcQ6ZJMWrjWwSbKlTZIkqQAc0yZJklQAhjZJkqQCMLRJkiQVgKFNkiSpAAxtkiRJBfB/+IiVT6KnMoAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 从数据得到参考周期图\n",
    "\n",
    "data = pd.read_csv(\"0.3_10_final.csv\")  \n",
    "\n",
    "dt=100\n",
    "counts_data = data['RATE']\n",
    "N=len(counts_data)\n",
    "pnum = np.arange(len(counts_data))\n",
    "t = [i*dt for i in pnum]\n",
    "\n",
    "'''\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(t,counts_data,'b')\n",
    "plt.xlabel(\"TIME(s)\")\n",
    "plt.ylabel(\"RATE(coount/s)\")\n",
    "plt.title(\"lightcurve\")\n",
    "plt.show()    \n",
    "'''\n",
    "\n",
    "nf = N/2\n",
    "df = 1/(dt*N)\n",
    "F_num = np.arange(1,nf)\n",
    "F = [i*df for i in F_num]\n",
    "mean_x = np.mean(counts_data)\n",
    "dft = fft(counts_data)\n",
    "dft1= dft[1:int(nf)+1]\n",
    "per_data = (abs(dft1)**2)*2*dt/((mean_x**2)*N)\n",
    "p_times_f_data = np.multiply(np.array(F1),np.array(per_data))\n",
    "\n",
    "\n",
    "\n",
    "# 数据分bin\n",
    "per_data_binned = databin_20(per_data)\n",
    "per_data_b_std = databin_20_std(per_data)\n",
    "'''p_times_f_data_b = np.multiply(np.array(F1_binned),np.array(per_data_binned))'''\n",
    "p_times_f_data_b = np.array(F1_binned)+np.array(per_data_binned)\n",
    "\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.scatter(F1_binned,p_times_f_data_b,color=\"r\",linewidth=1) \n",
    "plt.errorbar(F1_binned,p_times_f_data_b,yerr=per_data_b_std, fmt='.r')\n",
    "plt.plot(F1_log,PTF_ps_log,color=\"b\",linewidth=1) \n",
    "plt.xlabel(\"log frequency\")\n",
    "plt.ylabel(\"log power*frequency\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.921970124922109\n"
     ]
    }
   ],
   "source": [
    "# chi2 =（每一点500平均值 - 每一点数据得到）/ 每一点的500标准差  求和\n",
    "\n",
    "chi2=0\n",
    "for i in range(len(ADM_binned_ps)):\n",
    "    chi2 = chi2+(((ADM_binned_ps[i]-per_data_binned[i])/per_data_b_std[i])**2)\n",
    "\n",
    "print (chi2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2_of_model(A):\n",
    "    data = pd.read_csv(\"0.3_10_final.csv\")  \n",
    "    counts_data = data['RATE']\n",
    "    dt=100\n",
    "    mean_x_data = np.mean(counts_data)\n",
    "    N=len(counts_data)\n",
    "    N_randomlc=len(counts_data)\n",
    "    \n",
    "    omega = []\n",
    "    POW = []\n",
    "    DFT = []\n",
    "    fr = []\n",
    "    fi = []\n",
    "    f1 = []\n",
    "    f2 = []\n",
    "    p = []\n",
    "    f = []\n",
    "    \n",
    "    counts_list=[None for i in range(500)]\n",
    "    \n",
    "    \n",
    "    # 500条光变曲线\n",
    "    for a in range(500): \n",
    "        f_b=1.7E-4\n",
    "        alpha_H=3.8\n",
    "        alpha_L=1.0\n",
    "        for j in range(1,int(N_randomlc)+1):\n",
    "            omega.append(j/(N_randomlc*dt))\n",
    "            POW.append(((omega[-1]**(-alpha_L))/(1+(omega[-1]/f_b)**(alpha_H-alpha_L)))*A)\n",
    "            DFT.append(complex(np.sqrt(POW[-1]),np.sqrt(POW[-1])))\n",
    "            s1=np.random.normal(loc=0.0, scale=1.0, size=None)\n",
    "            s2=np.random.normal(loc=0.0, scale=1.0, size=None)\n",
    "            fr.append((DFT[-1].real)*s1)\n",
    "            fi.append((DFT[-1].imag)*s2)\n",
    "            f1.append(complex(fr[-1],fi[-1]))\n",
    "        counts = ifft(f1)\n",
    "        counts_list[a]=counts\n",
    "    \n",
    "    \n",
    "    #周期图\n",
    "    per_list=[None for i in range(500)]\n",
    "    for a in range(500):\n",
    "        nf = N/2 \n",
    "        df = 1/(dt*N)\n",
    "        F_a = np.arange(1,nf+1)\n",
    "        F = [i*df for i in F_a]\n",
    "        F1 = F[0:int(nf)]\n",
    "        mean_x = np.mean(counts_list[a])\n",
    "        x_new  = [i-mean_x for i in counts_list[a]]\n",
    "        dft   = fft(counts_list[a])\n",
    "        dft_1 = dft[1:int(nf)+1]\n",
    "        per = (abs(dft_1)**2)\n",
    "        per_list[a] = per\n",
    "    \n",
    "    \n",
    "    # 500条周期图取 mean,std\n",
    "    per_everypoint_list=[]\n",
    "    for i in range(int(nf)):\n",
    "        per_everypoint_list.append([])\n",
    "    ADM=[]\n",
    "    ADM_std=[]\n",
    "    for m in range(int(nf)):\n",
    "        for n in range(500):\n",
    "            per_everypoint_list[m].append(per_list[n][m])\n",
    "        \n",
    "    for m in range(int(nf)):\n",
    "        ADM.append(np.mean(per_everypoint_list[m]))\n",
    "        ADM_std.append(np.std(per_everypoint_list[m]))    \n",
    "        \n",
    "        \n",
    "    POW1=POW[0:int(nf)]\n",
    "    P_TIMES_F = np.multiply(np.array(F1),np.array(POW1))\n",
    "    F1_log = [math.log(i,10) for i in F1]\n",
    "    PTF_log = [math.log(i,10) for i in P_TIMES_F]\n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    #加上泊松项并分bin\n",
    "    ADM_with_poisson = [i+2/mean_x_data for i in ADM]\n",
    "    ADM_binned_ps = databin_20(ADM_with_poisson)\n",
    "    F1_binned = databin_20(F1)\n",
    "    ADM_times_f_b_ps = np.array(F1_binned)+np.array(ADM_binned_ps)\n",
    "    \n",
    "    \n",
    "    # 从数据得到参考周期图\n",
    "    pnum = np.arange(len(counts_data))\n",
    "    t = [i*dt for i in pnum]\n",
    "    \n",
    "    nf = N/2\n",
    "    df = 1/(dt*N)\n",
    "    F_num = np.arange(1,nf)\n",
    "    F = [i*df for i in F_num]\n",
    "    mean_x = np.mean(counts_data)\n",
    "    dft = fft(counts_data)\n",
    "    dft1= dft[1:int(nf)+1]\n",
    "    per_data = (abs(dft1)**2)*2*dt/((mean_x**2)*N)\n",
    "    \n",
    "    \n",
    "    # 数据分bin\n",
    "    per_data_binned = databin_20(per_data)\n",
    "    per_data_b_std = databin_20_std(per_data)\n",
    "    p_times_f_data_b = np.array(F1_binned)+np.array(per_data_binned)\n",
    "    \n",
    "    \n",
    "    chi2=0\n",
    "    for i in range(len(ADM_binned_ps)):\n",
    "        chi2 = chi2+(((ADM_binned_ps[i]-per_data_binned[i])/per_data_b_std[i])**2)    \n",
    "    \n",
    "    print('chi2 of this time = ',chi2)\n",
    "    \n",
    "    return (chi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 of this time =  2.7571222309627315\n",
      "2.7571222309627315\n"
     ]
    }
   ],
   "source": [
    "print(chi2_of_model(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 of this time =  2.3140334785244\n",
      "2.3140334785244\n"
     ]
    }
   ],
   "source": [
    "print(chi2_of_model(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 of this time =  2.762232826234981\n",
      "2.762232826234981\n"
     ]
    }
   ],
   "source": [
    "print(chi2_of_model(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 of this time =  2.323397419993069\n",
      "2.323397419993069\n"
     ]
    }
   ],
   "source": [
    "print(chi2_of_model(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "chi2 of this time =  2.9447632170961366\n",
      "2.9447632170961366\n"
     ]
    }
   ],
   "source": [
    "print(chi2_of_model(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.266"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.723\n",
    "2.633\n",
    "2.722\n",
    "2.356\n",
    "2.079\n",
    "2.266\n",
    "2.921\n",
    "2.757\n",
    "2.314\n",
    "2.762\n",
    "2.323\n",
    "2.944"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
